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Abstract 

Presented  is  a  multiscale  methodology  enabling  description  of  the  fundamental  mechanical 
behavior  of  crystalline  materials  at  the  length  scale  of  a  macroscopic  continuum  (i.e.,  millimeter 
resolution)  given  a  characterization  of  discrete  atomic  interactions  at  the  nanoscale  (i.e.,  angstrom 
resolution).  Asymptotic  homogenization  methods  permit  the  calculation  of  effective  mechanical 
properties  (e.g.  strain  energy,  stress,  and  stiffness)  of  a  representative  crystalline  volume  element 
containing  statistically  periodic  defect  structures.  From  a  numerical  standpoint,  the  theoretical- 
computational  method  postulated  and  implemented  here  in  the  context  of  lattice  statics  enables 
prediction  of  minimum  energy  configurations  of  imperfect  atomic-scale  crystals  deformed  to  finite 
strain  levels.  Numerical  simulations  demonstrate  the  utility  of  our  framework  for  the  particular  case 
of  body-centered-cubic  tungsten.  Elastic  stiffness  and  energetic  properties  of  periodic  unit  cells 
containing  vacancies,  screw  dislocations,  and  low-angle  twist  boundaries  are  computed.  Nonlinear 
aspects  of  elastic  behavior  in  the  context  of  plastic  flow  are  then  modeled  from  the  perspective  of 
atomistic-to-continuum  homogenization,  following  the  introduction  of  a  minimal  set  of  kinetic 
assumptions  required  to  account  for  the  propagation  of  dislocations  across  the  unit  cell  at  finite 
deformations. 
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1.  Introduction 

Atomistic-to-continuum  scale  transitioning  techniques  have  come  to  the  fore  in  recent 
years,  accompanying  and  supporting  widespread  interest  in  engineering  of  nanoscopic 
structures  (Rudd  and  Broughton,  2000;  Ghoniem  et  ah,  2003;  Fish  and  Chen,  2004). 
Although  marked  increases  in  computing  power  over  the  last  two  decades  have  facilitated 
large-scale  atomistic  calculations  (cf.  Seppala  et  ah,  2004),  spatial  and  temporal  scale- 
transitioning  schemes  are  presently  in  high  demand  since,  at  the  current  rate  of  increase  in 
computational  performance,  many  years  still  remain  before  computers  will  attain  the 
capability  to  simulate  the  mechanical  response  (and/or  thermal,  electrical,  or  magnetic 
response)  of  enough  atoms,  over  a  time  scale  of  relevant  duration,  to  realistically  address 
design  of  components  or  micro  structures  possessing  spatial  dimensions  on  the  order  of 
micrometers  or  larger. 

Existing  atomistic-to-continuum  modeling  methods  may  be  loosely  grouped  into  two 
categories:  simultaneous  or  coupled  approaches  and  sequential  or  decoupled  approaches. 
Simultaneous  methods  involve  concurrent  solution  of  governing  equations  at  multiple 
length  scales  (e.g.  quantum  mechanics,  molecular  dynamics,  continuum  mechanics),  with 
different  regions  of  the  spatial  domain  addressed  by  governing  equations  corresponding  to 
the  associated  scale  of  resolution  of  that  region.  Various  algorithms,  often  quite 
formidable  with  regards  to  numerical  implementation,  are  then  used  to  bridge  or 
“handshake”  neighboring  regions  of  differing  resolution,  for  example  regions  considering 
molecular  dynamics  and  continuum  mechanics.  Recent  examples  of  coupled  simultaneous 
methods  include,  for  example,  the  quasi-continuum  theory  (Tadmor  et  ah,  1996;  Shenoy 
et  ah,  1999),  the  three-scale  bridging  technique  of  Broughton  et  al.  (1999),  and  the  dynamic 
bridging  scale  method  of  Park  et  al.  (2005).  Such  methods  are  efficient  for  solving  problems 
wherein  the  bulk  of  the  material  remains  relatively  homogeneous  (and  amenable  to  a 
coarse-scale  solution),  with  the  hne  scale  physics  confined  to  a  small  volume  of  interest, 
such  as  at  a  crack  tip,  indentation  region,  or  grain  or  phase  boundary. 

In  contrast,  sequential  decoupled  methods  do  not  involve  direct  matching  of  boundary 
conditions  at  multiple  scales  of  resolution.  Instead,  fine  scale  calculations  are  conducted, 
and  then  the  pertinent  information  from  the  fine  scale  solutions  is  transferred  to  enhance 
the  representation  of  the  response  modeled  at  the  coarse  scale.  For  example,  atomistic- 
scale  simulations  of  the  motion  of  a  representative  number  of  atoms  may  be  used  to  extract 
bulk  elastic  and/or  plastic  properties  (Hao  et  al.,  2004)  or  interfacial  properties  (Spearot 
et  al.,  2005)  for  subsequent  implementation  in  continuum  mechanics-based  models. 
Methods  of  this  nature  include  approaches  based  on  simultaneous  conservation  of  mass, 
momentum,  and  energy  at  disparate  resolutions  (Zhou  and  McDowell,  2002)  as  well  as 
those  predicated  on  decomposition  of  the  displacement  vector  into  local  and  global 
components  (Hughes  et  al.,  1998;  Hao  et  al.,  2004). 

The  asymptotic  homogenization  technique  forwarded  in  the  present  paper  (see  also 
Chung  and  Namburu,  2003;  Chung,  2004)  falls  into  the  (latter)  category  of  spatially 
decoupled  methods.  In  this  approach,  discrete  calculations  are  conducted  at  the  atomistic 
level,  with  each  characteristic  volume  element  (i.e.,  unit  cell)  of  atoms  subjected  to  periodic 
boundary  conditions.  Asymptotic  homogenization  methods  (Bensoussan  et  al.,  1978; 
Sanchez-Palencia,  1980)  are  concurrently  employed  to  deduce  the  macroscopic  tangent 
stiffness  associated  with  the  mechanical  response  of  the  aggregate.  The  Cauchy-Born 
approximation  (cf.  Ericksen,  1984)  is  invoked  for  imposition  of  the  bulk  continuum 
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deformation,  with  the  fine  scale  deformation  of  the  atoms  identified  with  the  inner 
displacements  in  the  asymptotic  approximation.  In  this  way,  the  fine  scale  deformation  (of 
the  atoms)  near  defect  cores  is  corrected  with  the  inner  displacements  from  the 
homogenization  solution,  thereby  resulting  in  a  lower-energy  state  for  the  system  than 
would  be  achieved  via  imposition  of  a  uniform  deformation  gradient  over  all  atoms 
comprising  the  unit  cell.  The  present  approach  is  ideal  for  addressing  the  response  of 
microstructures  containing  spatially  periodically  distributed  defects  because  only  one  or  a 
few  defects  need  be  simulated  explicitly  at  the  atomistic  level  within  the  context  of  the 
periodicity  assumption  invoked  in  our  homogenization  scheme.  However,  due  to  this  very 
same  periodicity  assumption,  the  method  suffers  in  the  sense  that  isolated  (i.e., 
nonrepeating)  defects  cannot  be  easily  modeled.  Our  scaling  concept  is  sketched  in  Fig.  1. 

The  atomistic-continuum  homogenization  technique  forwarded  here  is  characterized  by 
several  distinctive  features  in  comparison  to  alternative  multiscale  formulations.  First, 
asymptotic  homogenization  methods  embody  weak  convergence  in  the  two-scale 
continuum  sense  (Nguetseng,  1989;  Allaire,  1992),  following  from  the  transmission  of 
displacement  gradients  as  opposed  to  displacements  across  length  scales.  Secondly,  this 
displacement-free  embedding  also  obviates  the  need  in  our  approach  for  potentially  overly 
restrictive  boundary  conditions  between  atomistic  and  continuum  domains.  Lastly,  in 
contrast  to  sequential  modeling  schemes  in  which  material  information  is  evaluated  at  a 
finer  scale  and  then  passed  to  a  coarser  one,  the  asymptotic  homogenization  method 
employs  a  self-consistent  set  of  governing  equations  simultaneously  involving  the  fine-  and 
coarse-scale  descriptions. 

Previous  research  (Chung  and  Namburu,  2003;  Chung,  2004)  addressed  the  effects  of 
pre-existing  point  defects  on  the  effective  elastic  stiffness  and  strain  energy  density  of 
graphene.  One  focus  of  the  present  effort  centers  on  the  mechanical  properties  induced  by 
dislocations,  the  primary  carriers  of  plastic  deformation  in  crystalline  metals.  Owing  to  the 
notion  that  the  flux  of  dislocations  leaves  the  crystal  lattice  relatively  unperturbed  upon 
complete  glide  across  the  periodic  unit  cell,  we  appeal  to  the  notion  of  an  evolving  relaxed 
intermediate  or  natural  configuration  (Eckart,  1948;  Bilby  et  ah,  1957;  Teodosiu,  1969; 
Asaro,  1983).  In  this  approach,  finite  inelastic  (i.e.,  plastic)  deformation  is  accounted  for  by 


Fig.  1.  Multiscale  modeling  schematic. 
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the  evolution  of  this  relaxed  conhguration,  as  dictated  by  the  dislocation  flux  through  the 
volume  element.  The  instantaneous  thermoelastic  response  of  the  element  is  measured  with 
respect  to  the  natural  conhguration,  and  driving  forces  for  defect  propagation  are 
calculated  at  the  current  state  or  time  instant,  in  order  to  project  the  intermediate 
configuration  forward  temporally.  Along  similar  lines,  in  an  effort  to  connect  nonlinear 
continuum  elasticity  theory  with  the  results  of  discrete  atomistic  calculations.  Till  and 
Broughton  (2000)  invoked  the  idea  of  an  ‘obviated  reference  lattice’  characterized  by  a 
slowly  varying  time  dependence  of  the  relaxed  intermediate  state,  with  a  strain  threshold 
introduced  above  which  relaxation  from  the  elastically  loaded  state  takes  place.  Groma 
and  Pawley  (1993)  investigated  an  analogous  problem  via  discrete  dislocation  dynamics 
simulations  (as  opposed  to  atomistics),  considering  in  their  case  the  superposition  of  elastic 
stress  flelds  from  the  defect  ensemble  as  well  as  the  energetically  driven  generation  of  new 
defects. 

Mechanical  properties  of  interest  in  the  present  context  include  the  elastic  stiffness  and 
energy  of  crystalline  materials  in  the  presence  of  defects  contained  within  the  volume 
element.  We  And  that  our  decoupled  approach  of  simultaneously  modeling  two  length 
scales  readily  enables  parametric  variations  of  the  defect  density  via  the  prescription  of  the 
number  of  atoms  in  the  flne  scale  representation  relative  to  the  total  number  of  defects 
embedded  within  the  periodic  unit  cell.  Our  framework  is  implemented  numerically  and 
applied  first  to  study  the  nonlinear  elastic  response  of  body-centered  cubic  (BCC)  tungsten 
(W)  containing  periodically  distributed  vacancies,  screw  dislocations,  screw  dislocation 
dipoles,  and  low-angle  twist  boundaries  (also  described  here  using  the  disclination  concept 
of  Li  (1972)).  In  the  above  simulations,  unit  cells  are  deformed  in  uniaxial  stretch  to  2.5% 
elongation.  Our  computational  method  is  demonstrated,  via  direct  comparison  with 
conjugate  gradient-based  lattice  statics,  to  be  an  efficient  means  of  predicting  minimum 
energy  configurations  of  imperfect  atomic-scale  lattices  subjected  to  finite  deformations. 
Next,  results  of  larger-deformation  shear  simulations  (to  10%  applied  shear  strain)  are 
reported  for  the  case  of  screw  dislocation  glide  in  BCC  W,  in  which  the  evolving 
intermediate  configuration  is  updated  assuming  monotonic  single  slip  occurs  as  resolved 
shear  stresses  exceed  the  Peierls  threshold  (cf.  Hirth  and  Lothe,  1982).  The  computational 
procedure  involves  determination  of  a  correction  to  the  macroscopically  imposed 
deformation  of  the  atoms  in  the  vicinity  of  defects.  Additional  details  regarding  numerical 
implementation,  efficiency,  and  validation  can  be  found  in  a  companion  paper  (Chung  and 
Clayton,  2006). 

The  remainder  of  this  paper  is  organized  as  follows.  Section  2  features  derivations  of  the 
multiscale  homogenization  equations  applicable  to  nonlinear  elastic  problems.  Section  3 
introduces  the  concept  of  an  evolving  reference  configuration  to  address  finite  deformation 
plasticity  and  the  corresponding  inelastic  deformation  gradient,  as  well  as  transformation 
rules  for  the  multiscale  equations  in  the  presence  of  this  evolving  frame.  Additional 
assumptions  and  limitations  of  our  method  relating  atomistic  variables,  defect  kinetics, 
and  continuum  flelds  are  given  in  Section  4,  followed  by  a  description  of  the  numerical 
implementation  in  Section  5.  We  present  results  of  demonstrative  simulations  in  Section  6 
and  offer  perspective  on  the  utility  of  our  theory  and  results  in  Section  7.  Vectors  and 
tensors  are  written  in  boldface  type,  with  scalars  and  individual  components  of  vectors  and 
tensors  written  in  italic  font.  The  indicial  notation  is  frequently  employed  for  clarity  of 
presentation,  with  summation  implied  over  repeated  indices,  e.g.  A'Bj  —  A^B\+  A^Bi  + 
A^B},  when  i=  1,2,3. 
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2.  Two-scale  asymptotic  homogenization 


2.1.  Governing  equations 


Reference  and  current  configurations  of  a  continuous  body  are  introduced,  denoted  by 
Bq  and  B,  respectively.  Let  X  and  x  represent  coordinates  spanning  the  reference  and 
spatial  frames,  and  let  x“  =  t),  with  t  denoting  time.  The  deformation  gradient  or 

tangent  mapping  F  from  Bq  to  B  is  then  written 


F  = 


0X 


 0X“ 


(1) 


The  following  strain  measures  are  also  introduced,  where  = 
and  g^i,  —  0aX  •  d/jX  denote,  respectively,  metric  tensors  in  reference  and  spatial  coordinate 
systems: 


Cab  —  pAdabC^B^  2Eab  —  F^AQahF^.B  ~  Cab- 


(2) 


Let  r,  P,  and  S  denote  the  Cauchy  stress,  first  Piola-Kirchhoff  stress,  and  second 
Piola-Kirchhoff  stress,  respectively,  related  by 


^ah  ^  J-l  palpi’ A  ^ 


(3) 


Localized  forms  of  the  balances  of  linear  and  angular  momentum  are  written  as  follows 
in  the  reference  frame,  assuming  quasi-static  conditions: 


nciA 

^\A 


B“^0, 


pa  pbA  _  paAph 


(4) 


where  the  vertical  bar  denotes  covariant  differentiation  and  B  is  the  body  force  vector  per 
unit  reference  volume.  From  Eq.  (3)  and  the  second  of  Eq.  (4),  we  have  Z"*  =  and 
where  parentheses  indicate  symmetrization,  i.e.,  —  A‘‘^  +  A^“  for 

arbitrary  second-rank  tensor  A.  As  will  be  discussed  later,  neglect  of  inertia  in  Eq.  (4) 
restricts  our  present  approach  to  time-independent  applications  at  the  fine  scale  (i.e.,  static 
problems).  Multiplying  the  first  of  Eq.  (4)  by  virtual  displacement  and  integrating  over 
reference  volume  V,  we  arrive  at  the  following  virtual  work  principle  (cf.  Marsden  and 
Hughes,  1983): 


[  P“^g,i,i3u)';j,dV^  [  T‘'g,,8u^dA+  [  Bl'g,M  dV,  (5) 

Jv  JdV  Jv 

with  the  traction  per  unit  reference  area  A  given  by  7“  —  P“^Nb,  where  are  covariant 
components  of  the  unit  normal  vector  to  external  boundary  0F.  Let  us  assume  the 
existence  of  a  free-energy  potential  P  per  unit  reference  volume  on  Bq,  with  the  stress 
tensor  satisfying  the  following  hyperelastic  relationships: 


SAB 


2 

^Cab 


QEab 


(6) 


For  the  particular  case  of  first-order  hyperelasticity,  Eq.  (6)  becomes 


SAB 


dEAB^EcD 


r'  _  (pABCD  rp 

^CD  — 


(7) 
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where  —  £(dB)(CD)  fourth-rank  elastic  modulus  tensor  in  the  reference  frame. 

Substituting  Eq.  (6)  into  Eq.  (5)  and  appealing  to  Eq.  (3),  we  then  have 


2.2.  Asymptotic  homogenization 


Our  link  between  atomistic  and  continuum  resolutions  is  achieved  via  the  asymptotic 
homogenization  technique  (Chung  and  Namburu,  2003;  Chung,  2004).  Eet  fine  and  coarse 
length  scales  (i.e.,  the  microscale  and  macroscale)  be  spanned  by  coordinates  y"  =  y“{Y^,  t) 
and  —  x‘‘{X"^,t),  respectively.  Notice  that  in  the  present  scheme,  both  scales  are 
parameterized  by  the  same  temporal  variable  t.  Multiscale  coordinates  are  related  by 

x"  =  e/,  X'^^eY^,  (9) 

where  a  is  a  small  scalar  assumed  constant  throughout  the  time  history  of  deformation.  We 
introduce  coarse-  and  hne  scale  displacements  u  and  v,  respectively,  restricted  below  to 
coincident  Cartesian  coordinate  systems  in  the  reference  and  spatial  frames: 

u“  =  x"  -  i;“  =  /  -  8“^  Y^,  (10) 

with  the  Cartesian  shifter  8“^  —  1  for  a  —  A  and  <5“^  =  0  for  Deformation  gradients 

then  follow  as 


F‘'  = 


8m" 


+  ^“a,  r.A^ 


0>>“ 


+  d 


(11) 


dXA  '  dY^  ' 

Next  we  assume  an  additive  decomposition  of  displacements  at  the  coarse  scale  (Takano 
et  ah,  2000), 


.  .a  TlU  I  ~a  -r-A  I 

U  =  U  +  W  =  U  +£1), 


(12) 


where  m“  represents  the  displacement  that  would  exist  in  a  microscopically  homogeneous 
medium  and  m“  is  the  perturbation  in  displacement  due  to  fine  scale  heterogeneity,  with 
corresponding  fine  scale  representation  v“.  The  corresponding  microscopic  decomposition 
is 


1  ^  ,a  -^a  I  ~a  (  rpa  ca  \  vA  i  c-.a 
E  U  —  V  +  V  —  ^  —  8^)  Y  +v, 


(13) 


with  the  microscopic  displacement  arising  from  the  projection  to  the  fine  scale  of  the 
macroscopic  deformation  gradient  F“^.  Differentiating  u  of  Eq.  (12)  with  respect  to  X'^  gives 


0 


(m"  -t-  eC“)  = 


0m“  ^  0)5“ 


(14) 


where  we  have  appealed  to  the  second  of  Eq.  (9).  Erom  Eq.  (11),  we  then  arrive  at  the 
decomposition 


0m“  dv 


=  Kf‘ 


A’ 


(15) 


where  F^'-  is  the  deformation  gradient  under  microscopically  homogeneous  conditions  and 
F’h  =  F  depends  upon  the  gradient  of  v  and  thus  accounts  for  heterogeneity  due  to 
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defect  fields.  Barred  Greek  indices  denote  a  third  configuration  B  for  the  material  associated 
with  the  multiplicative  decomposition  in  the  last  of  Eq.  (15). 

Note  that  the  left-hand  side  of  Eq.  (8)  can  be  written  as  follows  in  Cartesian  coordinates: 

d>P  0((5m") 


f  d'P 

'V  OCab 


dV, 


IvdF^B 

and  that  the  total  displacement  variation  3iF  can  be  expressed,  from  Eq.  (12),  as 
Sif  —  dtf  +  sSv". 

Snbstitnting  Eqs.  (16)  and  (17)  into  Eq.  (8)  gives 
/•  05'  8 

lvd± 


(16) 


(17) 


■0E“  0Z^ 


(du"  +  eSv^)  dV 


=  f  r'g^i,i5u^  +  s5v>’)dA+  f  B‘'gJSJl^  +  s3v>’)dV, 
JdV  Jv 


(18) 


IdV 

which  is  then  volume-averaged  over  micro-domain  Y  to  yield 

1 

Y 


f  f  1 

fd(3u‘‘)  , 

0(^i;")\ 

JYJvdF^j,' 

[dx^  ' 

07^  ) 

dVdY 


Eq.  (19)  is  satisfied  in  the  asymptotic  limit  £->0  only  if  the  following  conditions  hold: 


1 


dF  0((5m“) 


YJYJvdF^B  dX^ 


dFdE  = 


/0F 


T‘‘g^h^u'’  dA 


1 


YJYJvdF^B  07 


"^dFd7  =  o  (Vdro. 


B“g^h3u‘’dV  (ySu^),  (20) 

(21) 


Notice  that  coarse-scale  Eq.  (20)  and  fine  scale  Eq.  (21)  are  coupled  through  the 
constitutive  dependency  W  —  'f'(F(0u/0X,  0v/0Y)). 


3.  Elomogenization  under  an  evolving  intermediate  frame 

3.1.  Multiplicative  kinematics 

Assume  now  that  we  are  dealing  with  an  elastic-plastic  material  in  the  usual  continuum 
sense  wherein  the  deformation  gradient  F  is  decomposed  multiplicatively  into  a  lattice  part, 
F^,  and  a  plastic  part,  F**  (Bilby  et  ah,  1957;  Teodosin,  1969): 

F  =  F^F**,  F  =  F^M^Y  (22) 

Implicit  in  Eq.  (22)  is  the  existence  of  intermediate  configuration  B,  with  F**  the  tangent 
mapping  between  B^  and  B,  and  F^  the  tangent  mapping  from  B  to  B.  The  plastic  map  F** 
reflects  contributions  to  the  total  deformation  that  leave  the  lattice  unperturbed,  such  as 
the  flux  of  mobile  dislocations.  The  lattice  map  F^  inclndes  all  other  kinematic 
contribntions,  including  rigid  body  motion  of  the  lattice,  recoverable  elastic  stretch  fields 
associated  with  the  applied  loads  on  the  body,  and  residual  elastic  stretch  fields  attributed 
to  microscopic  heterogeneity  (e.g.,  defects  contained  within  the  local  volume  element  or 
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periodic  unit  cell).  Thus  B  is  identified  with  a  relaxed  natural  configuration  of  finite 
plasticity  theory  (Eckart,  1948). 

3.2.  Governing  equations 


In  an  elastic-plastic  body,  it  is  typically  assumed  that  the  free  energy  depends  only  upon 
that  part  of  the  deformation  that  affects  the  lattice,  in  other  words, 

ip  =  T(F,  pP,  I,  B)  =  •T(F^  1, 0),  (23) 

where  P  is  measured  per  unit  intermediate  volume  on  B,  is  the  Jacobian  deter¬ 
minant  of  F*”,  and  where  we  have  also  included  the  dependence  of  the  free  energy  on 
absolute  temperature  0.  Additionally,  |  is  a  vector  of  internal  state  variables  accounting 
for  deviations  in  stored  energy  from  that  of  a  perfect  lattice.  From  Eqs.  (22)  and  (23) 
and  assuming  an  instantaneous  hyperelastic  response,  the  first  Piola-Kirchhoff 
stress  satisfies 


P 


.A 

a 


_  jp-i  8^  pP-lA 
.a 


The  balance  of  linear  momentum  in  Eq.  (4)  can  then  be  re-written  as  follows: 


(24) 


dP  \ 


+  B“  =  0. 

\A 


Eet  us  denote  the  intermediate  second  Piola-Kirchhoff  stress  by 


dP 


dF 


Lh 


pL-lfS 


dP  _  dP 

dC^ls  dE^fi  ’ 


where 


(25) 


(26) 


Cap  =  F^gahF^p\  '^Eap  =  Q/j  -  g^p,  (27) 

with  g  the  metric  tensor  on  configuration  B.  Note  that  in  Eq.  (26),  we  have  made  the  usual 
assumption  that  the  dependency  of  P  upon  F^  can  be  replaced  with  a  dependency  upon  C 
or  E,  in  order  to  satisfy  standard  frame  indifference  arguments  (cf.  Marsden  and  Hughes, 
1983).  Analogous  to  Eq.  (7),  for  the  particular  case  of  first-order  hyperelasticity,  Eq.  (26) 
becomes 


dEapdE^s 


E,s  = 


(28) 


with  —  £GP)(xS)  fourth-rank  elastic  modulus  referred  to  B.  Eq.  (5)  can  then  be 
written  as  follows  for  a  hyperelastic-plastic  body,  upon  appealing  to  Eq.  (25) 


/  /^-‘F^,-'^^(<5n)t4dF=  [  T‘'g,,5u^dA+  [  dV. 

JV  OP  a.  70F  Jv 


(29) 
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3.3.  Asymptotic  homogenization 


Expressions  (9)-(15)  and  (17)  apply  for  the  elastic-plastic  case.  Substituting  Eq.  (17)  into 
Eq.  (29)  and  using  Cartesian  coordinates  we  have 


J 


/._!  0  ^  e5v“)dV 


=  f  r'g^hi5u‘’ +  £5v>’)dA+  f  B“g^^iSu’’  +  sdv^)  dV,  (30) 

JdV  Jv 

leading  to  the  following  coupled  multiscale  equations  in  the  limit  as  £  ^  0: 


jP-lpP-lB  ^y^Y: 

dF^  dX^ 


IdV 


T^g^h^WdAF  /  B“g„i,8W  dV  (V  hM"), 
Jv 

(31) 


1 

Y 


jP-l 


dF  diSv‘‘) 

QpLa  QyB 


dFdE=0 


(V(5C“). 


(32) 


Notice  that  generally,  and  F*”  of  Eq.  (22)  and  F  and  F  of  Eq.  (15)  are  four  distinct 
deformation  mappings,  with  corresponding  conhgurations  of  the  body  depicted  in 
Fig.  2(a).  The  former  two  denote  lattice  and  plastic  deformations;  the  latter  two  represent 
micro-homogeneous  and  micro-heterogeneous  deformations,  respectively.  Under  certain 
circumstances,  however,  some  of  these  may  reduce  to  the  identity  map.  Possible  scenarios 
are  summarized  below: 


Elastoplasticity  with  defects  :  F  =  F^F**  =  FF, 

Homogeneous  elastoplasticity  :  F  =  F^F**  =  F, 

Elasticity  with  defects  :  F  =  F^  =  FF, 

Homogeneous  elasticity  :  F  =  F^  =  F.  (33) 

In  Eq.  (33),  “elastoplasticity  with  defects”  is  the  most  general  case,  for  example  a 
crystalline  volume  element  that  has  sustained  dislocation  flux  during  its  deformation 
history  and  also  contains  dislocations,  possibly  immobile.  By  “homogeneous  elastoplas¬ 
ticity”  we  mean  a  situation  in  which  dislocations  may  have  traversed  the  volume  element. 


Fig.  2.  Deformation  mappings  and  configurations  (a),  and  purely  plastic  deformation  (b). 
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but  no  defects  or  commensurate  residual  stress  fields  are  contained  within.  By  “elasticity 
with  defects”,  we  refer  to  the  case  in  which  no  defects  have  traversed  the  element  (i.e.,  null 
dislocation  flux  over  the  deformation  history),  but  when  defects  are  contained  within, 
leading  to  microscopic  heterogeneity  and  internal  stress  fields.  This  was  the  situation 
studied  by  Chung  and  Namburu  (2003)  and  Chung  (2004),  where  the  representation  of 
point  defects  was  the  primary  concern.  “Homogeneous  elasticity”  refers  to  a  lattice  free  of 
defects  over  the  entire  time  history  of  deformation. 

4.  Multiscale  formulation 

Here  we  present  kinematic,  thermodynamic,  and  kinetic  assumptions  needed  to  relate 
atomistic  and  continuum  fields  and  incrementally  update  intermediate  configurations  B 
and  B. 


4.1.  Discrete  and  continuum  descriptions  of  lattice  quantities 


Assume  that  in  reference  configuration  Bq,  the  representative  volume  (i.e.,  unit  cell)  for 
homogenization  consists  of  atoms  arranged  in  a  lattice,  perhaps  imperfect  due  to  the 
presence  of  defects.  Furthermore,  assume  that  in  configuration  B,  the  same  mass  and 
number  of  atoms  exist  in  this  representative  volume,  for  example  upon  propagation  of 
dislocations  across  the  unit  cell.  Individual  atomic  positions  in  Bq  and  B  may  not  coincide, 
even  though  the  lattice  may  look  the  same  to  an  external  observer  in  each  of  these  two 
configurations.  The  position  vector  for  an  arbitrary  atom  m  in  evolving  configuration  B  is 
denoted  by  Z(„) ,  where  angled  brackets  are  reserved  for  atomic  labels,  which  range  from  1 
to  I\l.  Spatial  positions  of  atoms  in  configuration  B  are  then  found  as  follows,  in 
Cartesian  coordinates: 


(34) 


with  a  displacement  vector  between  intermediate  and  current  states  for  atom  m.  Let 
Rij\k)  and  rijyk)  denote  vectors  separating  atoms  a  and  b  in  respective  configurations  B  and 
B,  i.e.. 


R 


m 


—  ^{k)  —  Z 


W’ 


(35) 


tt{]\k)  —  Z(/c)  Zy).  (36) 

Making  contact  with  the  homogenization  theory  of  Sections  2  and  3,  we  next  assume 

+  (37) 

with  summation  implied  over  repeated  atomic  indices.  The  first  term  in  Eq.  (37),  , 

accounts  for  the  uniform  projection  over  each  periodic  cell  of  the  macroscopic  lattice 
deformation  field  to  the  fine  scale  (i.e.,  the  Cauchy-Born  approximation),  and  is  the 
discrete  atomistic  analog  of  the  perturbation  in  displacement  due  to  microscopic 
heterogeneity  given  in  previously  in  Eq.  (13),  written  here  for  atom  m.  Notice  that 
of  Eq.  (37)  is  technically  the  perturbation  in  displacement  from  the  intermediate  state, 
whereas  v"  of  Eq.  (13)  is  a  displacement  perturbation  from  the  initial  reference  state. 
However,  if  the  plastic  deformation  mapping  F**  is  assumed  to  leave  the  lattice  properties 
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unaltered,  then  these  displacement  perturbations  may  be  considered,  for  our  purposes,  to 
be  equivalent.  Spatial  separation  vectors  in  Eq.  (36)  then  become 

^{f\k)  =  +  ^(k)  -  y(j')  =  (38) 

with  iijyk)  —  V((t)  —  V(y)  accounting  for  deviations  from  the  Cauchy-Born  rule. 

Reference  state  Bo  is  related  to  the  intermediate  state  B  as  follows.  The  mapping  F** 
represents  the  cnmulative  deformation  of  the  material,  but  not  the  lattice,  due  to 
dislocation  glide  through  the  unit  cell.  For  a  hxed  control  volume,  once  atoms  have 
convected  through  the  lattice  due  to  F**,  Bq  and  B  appear  identical  in  the  fine  scale  domain 
if  no  dislocations  are  created  or  destroyed  within  the  volume  element,  yet  net  deformation 
of  the  material  will  have  taken  place  at  the  coarse  scale.  The  same  concept  is  used  in 
classical  crystal  plasticity  theory,  e.g.  Asaro  (1983),  wherein  the  lattice  director  vectors  are 
assumed  to  be  unaffected  by  F**.  For  a  purely  plastic  process  with  F^  =  1,  new  atoms 
would  enter  the  control  volume  (i.e.,  the  atomic  unit  cell)  in  the  identical  locations  as  the 
old,  to  replace  those  that  exited  due  to  plastic  flow,  and  the  mass  of  the  system  would  be 
conserved  (the  volume  would  also  be  conserved  for  processes  in  which  —  1).  The 
concept  is  illustrated  in  Fig.  2(b),  in  which  atoms  that  exit  the  unit  cell  are  denoted  by  filled 
circles,  and  those  that  enter  the  cell  by  open  circles,  with  the  underlying  material 
undergoing  pnre  plastic  shear.  Though  not  illnstrated  explicitly  here,  immobile 
dislocations  or  other  stationary  defects  are  admitted  within  the  volume  in  both 
configurations,  as  the  influence  of  such  defects  on  the  stored  energy  and  mechanical 
response  is  of  primary  interest  in  the  applications  that  follow. 

Flenceforth  we  assume  a  free-energy  potential  depending  only  upon  the  relative 
positions  of  all  atoms  (i.e.,  lattice  statics): 

^  =  ^'(r(i\2),r(i\3),i‘(2\3)>  ■  ■  ■  r(N-i\N))-  (39) 

Analogous  to  Eq.  (39)  is  the  isothermal  representation  of  continuum  expression  (23): 

^'^nFM),  (40) 

where  |  accounts  for  effects  on  stored  energy  due  to  deviations  from  homogeneous 
elasticity  at  the  fine  scale.  Integrating  Eq.  (32)  by  parts  and  applying  the  divergence 
theorem  over  volume  Y  with  oriented  surface  element  dA,  we  arrive  at 


1 

A 


f  / 

JdY  Jv 


jP-^pP-^B 


dW 

dF^' 


Sv^NBdVdA 


dFdT. 

(41) 


Focalizing  the  volume  integral  on  the  right-hand  side  of  Eq.  (41)  and  considering  all 
admissible  variations  5v,  the  microscopic  linear  momentum  balance  becomes 


^  f  rP-lpP-l  B  \ 

0T^V  ■“ 


dPl 

0T^ 


=  0 


(m  F) 


(42) 


as  the  area  integral  vanishes  since  we  may  choose  =  0  on  ^  =  0F.  Eq.  (42)  is  quite 
general  in  the  sense  that  no  assumption  is  made  on  the  order  of  the  incremental  elastic 
response;  for  example,  higher-order  elastic  constants  are  admitted  in  the  sense  of  material 
nonlinearity.  For  algorithmic  purposes,  however,  we  find  it  advantageous  to  assume  a 
quadratic  energy  dependency  along  the  lines  of  Eq.  (28),  an  appropriate  assumption  for 
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most  engineering  metals  undergoing  quasi-static  deformations.  The  overall  elastic 
stress-strain  response  may  still  be  nonlinear  in  our  formulation,  however,  as  we  allow 
the  tangent  elastic  stiffness  tensor  to  vary  with  elastic  deformation  F*",  and  the  formulation 
may  also  be  considered  nonlinear  in  the  geometric  sense,  as  hnite  deformations  are 
addressed.  We  thus  re-write  Eq.  (42)  as 


0 


(43) 


where  C^^,,  is  a  mixed-variant  effective  elastic  modulus  tensor.  Applying  our  previous 
assumptions  that  F**  does  not  affect  the  lattice  state  and  acts  uniformly  over  Y,  Eq.  (43) 
becomes 


jP-lpP-lB 


0r 


(44) 


Next,  from  the  deformation  gradient  expression  given  by  the  hrst  of  Eq.  (15),  we  have 


p^lB^^^ah 


0F"  V0T^' 


,77* 


du‘ 


^.A 


s’’  F^} 

■X 


^P-\A  _  pP-\B 


07^ 


ah 


07" 


■fP-IA 


(45) 


Our  immediate  goal  is  to  express  Eq.  (45)  in  terms  of  atomistic  and  macroscopic 
displacements.  Invoking  the  chain  rule  gives 


0  ^07-^  0  ^yP-\A  0 

0Z"^.)  “0Z“)07-^  “  07^^ 


(46) 


_ P—\ji 

where  the  linear  operator  F^j^^  measures  the  change  in  intermediate  atomic  position  of 
atom  j  with  respect  to  a  change  in  reference  coordinates.  Efsing  Eq.  (46)  in  the  left-hand 
side  of  Eq.  (45),  we  assume 


^P-\B^^ah 

■“  07^ 


,77* 


du‘ 


+  S' 


dX'^ 
f  du’’ 


.A 


s’’ 

-A 


^se/? 


^P-IA 


0c„;,/0M' 


0Z,) 


Px 


A  +  ^  A  -  ^"yF.A 


pP-\A 


177* 


V0z^ 


s’’  F^^- 

■X  -A 


7P-IA 


f- 


377* 


</■)«*  \q^A 


+  <5  .  -  KF 


'Px 


pP-\A 


(47) 


and  for  the  right-hand  side  of  Eq.  (45),  rigorous  for  small  atomic  perturbations 


yP-lB 


07^ 


^ah 


dv’’ 


jyP-\A 


0Z 


^ah  ~~!r 


(j) 


dZ‘ 


’’{k} 


^{jk)ah^{k)- 


(48) 


In  Eqs.  (47)  and  (48),  the  notation  denotes  the  transformation  steps  v”  ^  and 
^  latter  implying  that  acts  uniformly  over  all  atoms  k 

within  each  periodic  cell.  Using  Eqs.  (47)  and  (48),  the  fine  scale  equilibrium  equation  (42) 
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finally  becomes 

{j)ah  I 


+  .5' 


.A 


-A 


^P-IA 


=  je 


{ik)ahli(k)  ’ 


(49) 


where  and  the  Hessian  matrix  =  0^^'/0?“^;)0<?^)-  In 

subsequent  calculations,  Eq.  (49)  is  solved  for  the  (atomic)  inner  displacements 
Notice  that  for  the  particular  case  of  an  elastic  (as  opposed  to  elastic-plastic)  material  with 
defects,  we  have  the  conditions  F  =  F^,  F**  =  1,  ?'(F^,|)  ^  'F(F,  |), 

^{jk)ah  —  nnd  Eq.  (49)  becomes  equivalent  to  the  original  proposition 

of  Chung  and  Namburu  (2003),  i.e., 

/  0j^/>  \ 

—  ^(jk)ahv\k)-  (50) 


Now  we  reconsider  coarse-scale  Eq.  (31),  assuming  that  v  is  known  from  solution  of 
Eq.  (49).  The  left-hand  side  of  Eq.  (31)  can  be  written  as 
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dFdT,  (51) 


where  fine  scale  perturbations  and  atomic  positions  affect  F^^  and  through  Eqs.  (15) 
and  (39),  respectively,  and  hence  influence  the  macroscopic  response.  Eollowing 
Eqs.  (46)-(49),  we  may  write 
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(52) 


such  that  Eq.  (31)  becomes 


ha  \^g^.4  -A  .X  .A  j  0^5 


=  T‘'g,M  dA  +  dV  +  j  dFdT.  (53) 


where  and  Later,  Eq.  (53)  will  be 

converted  to  a  finite-element  equation  applied  over  spatially  discretized  domain  F. 

As  is  clear  from  relations  (47)  and  (48),  the  theory  effectively  equates  f"  and  ^“(kP  implying 
two  main  assumptions.  Firstly,  neglecting  the  change  in  configurations,  the  atomic  discrete 
degrees-of-freedom  are  assumed  to  be  equivalent  to  the  fine  scale  continuum  analog  from 
classical  two-scale  homogenization  theory,  f®.  This  assertion  was  made  in  earlier  papers  by 
Chung  and  co-workers  (Chung  and  Namburu,  2003;  Chung,  2004),  and  was  further  modified 
for  efficient  computer  implementation  by  the  present  authors  (Chung  and  Clayton,  2006).  It 
is  a  fundamental  assumption  made  without  mathematical  proof.  Physically,  it  implies  that 
atomic  deviations  from  the  Cauchy-Born  rule  are  equivalent  to  fine  scale  perturbations  due 
to  microstructure  heterogeneity  in  the  sense  of  two-scale  continuum  homogenization. 
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Secondly,  u"  is  a  spatial  quantity  in  the  two-scale  continuum  theory  representing  local 
perturbations  from  the  reference  state  Bq,  whereas  is  a  spatial  quantity  representing 
atomic  perturbations  from  the  intermediate  configuration  B.  Since  F**  does  not  affect  the 
underlying  lattice  or  the  atomic  coordinates  as  discussed  above,  we  assume  that  the 
perturbations  from  the  reference  state  are  the  same  as  those  from  the  intermediate  state. 

The  manipulations  in  Eq.  (47)  effectively  embed  a  deformation  gradient  (i.e.,  F^)  into 
the  operator  Hence  the  rightmost  term  of  the  macroscopic  equation  (53)  implicitly 

involves  gradients  of  the  atomistic  variable  This  result  follows  from  the  two-scale 
mathematical  homogenization  approach  presently  used  and  embodies  the  features  of  weak 
convergence  usually  encountered  for  such  methods  (Nguetseng,  1989;  Allaire,  1992),  as 
discussed  earlier.  In  practice  and  in  contrast  to  many  existing  coupled  multiscale 
approaches,  the  characteristic  length  scales  in  X  (continuum)  and  Y  (atomistic)  are  not 
necessarily  identical,  and  therefore,  no  length  scale-preserving  displacement  boundary 
conditions  between  the  domains  are  required.  Only  the  correlated  inter-scale  gradients  are 
enforced  in  an  average  sense  over  each  domain.  We  also  note  that  Eqs.  (49)  and  (53)  are 
coupled  in  terms  of  solution  variables  v  and  u. 

4.2.  Inelastic  fields  and  defect  kinetics 

Depending  upon  the  fine  scale  physics  representation,  additional  kinetic  equations  may 
be  required  in  our  framework  in  order  to  advance  the  plastic  deformation  gradient  F**  and 
internal  variables  |.  In  the  continuum  description,  these  are  written  as 


f'’  =  f'’(f,fM,0), 

(54) 

|  =  |(f,fM,0), 

(55) 

assuming  that  instantaneous  rates  depend  upon  the  thermodynamic  state  described  by  the 
first  of  Eq.  (23).  When  a  finite-temperature  atomistic  description  is  invoked  at  the  fine 
scale,  Eq.  (54)  may  not  be  needed  a  priori,  as  the  plastic  velocity  gradient  can  conceivably 
be  expressed  in  terms  of  the  flux  of  mobile  dislocations.  One  can  formulate  an  approximate 

•  P 

expression  for  the  rate  of  the  mapping  F  in  terms  of  the  dislocation  velocity,  line  length, 
and  orientation  (Teodosiu,  1969).  A  similar  approach  has  been  used  by  Zbib  and  co¬ 
workers  in  the  context  of  discrete  dislocation  plasticity  (Zbib  and  de  la  Rubia,  2002). 

.  p 

Conceivably,  our  approach  could  be  modified  such  that  F  could  be  extracted  in  such  a 
manner  directly  from  the  atomistics,  though  it  is  a  nontrivial  task  to  precisely  define  the 
location  of  the  dislocation  line  in  the  atomistic  domain,  and  our  method  would  need  to  be 
extended  to  account  for  dynamics  in  order  to  track  the  dislocation  velocity,  as  no  time 
scale  enters  the  present  static  formulation.  One  cannot  construct  an  explicit  formula  for  the 
tangent  map  F**  itself  unless  the  dislocation  geometry  and  velocity  are  constant  over 
the  time  history  of  deformation.  In  the  present  implementation,  Eq.  (54)  is  needed  as  the 
physics  of  thermally  activated  dislocation  glide  are  not  adequately  addressed  by  our  static 
calculations  at  the  fine  scale.  Relation  (55)  controls,  for  example,  rates  of  defect  generation 
or  annihilation,  in  instances  when  such  kinetics  are  not  modeled  directly  at  the  fine  scale,  as 
is  again  the  case  in  the  present  zero-temperature  atomistic  implementation.  Simple, 
yet  specific,  forms  of  Eqs.  (54)  and  (55)  follow  in  the  discussion  of  numerical  results  in 
Section  6.3. 
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5.  Numerical  implementation 


As  the  material  response  is  nonlinear  in  the  presence  of  defects,  an  iterative  scheme  is 
employed  for  application  of  our  computational  method  to  deforming  crystals.  The  initial¬ 
boundary  value  problem  consists  of  the  macroscopic  domain  V,  discretized  here  into 
standard  hnite  elements.  The  integration  point  of  each  element  in  turn  represents  a 
microscopic  unit  cell  of  volume  Y,  consisting  of  N  atoms  subjected  to  periodic  hne  scale 
boundary  conditions.  Defects  such  as  dislocations  may  be  present  in  the  initial 
conhguration,  at  ?  =  0,  accompanying  the  standard  assumption  F^q  =  1.  Details  of  the 
iterative  algorithm  are  listed  below.  In  what  follows,  we  denote  by  At  the  time  increment 
over  which  a  numerical  integration  cycle  takes  place,  and  AO  and  Av  denote  increments  in 
macro-  and  micro-displacement  vectors  over  At. 


1.  Incrementally  increase  the  applied  load. 

(a)  Begin  element  loop,  for  each  element  level  integration  point: 

(i)  Update  the  coarse-scale  deformation  gradient  F. 

(ii)  Compute  the  lattice  deformation  gradient  F^  =  FF**"',  where  at  the  coarse 
scale,  F  and  F  are  assumed  effectively  equivalent. 

(hi)  Update  atomic  positions  via  Eq.  (37),  with  v  =  0  initially  and 

(iv)  Calculate  the  Hessian 

from  the  atomistic  potential  energy  E  (for  example,  see  later  Eqs.  (70)-(71)). 

(v)  Apply  periodic  boundary  conditions  and  solve  fine  scale  Eq.  (49)  for  Av: 


(J)ah  1 


fdiAu>^  i, 


s’’ 

.7  .A 


:h  'V 
(k)) 


(56) 


where  Au  is  known  at  present  following  the  update  of  F  of  step  (i)  above. 

(vi)  Repeat  steps  (iii),  (iv),  and  (v)  above  until  Av  ^  0  on  successive  counts. 

(vii)  Compute  /dF^^dF^I’  from  the  atomic  energy  E  (see  later 


Eq.  (72)). 


(viii)  Calculate  the  effective  elastic  stiffness  and  the 


quantity  =  J 


P-\ 


(Dha- 


(b)  Compute  and  assemble  the  finite  element  equations.  In  discretized  incremental 
form,  the  left-hand  side  of  (53)  becomes 
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(57) 


where  N  is  the  element  shape  function  and  Au  is  the  nodal  displacement 
increment.  Similarly,  the  right-hand  side  of  (53)  can  be  written 


AT‘‘g,MdA+  I  A5“g,,AM'’dF+- 
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IdV 
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0 

Jv  ”■  "'0X 

Combining  (57)  and  (58)  then  gives 
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(59) 


which  can  be  reduced  upon  taking  the  hrst  variation  to  the  hnite  element  equation 
K,.gAt/  =  Z,  +  Z,  +  D^-  (60) 

Solution  of  Eq.  (60)  involves 

(i)  Compute  —  y^,c  ^y\d(N^^) ldX^'\J^\d(Nj) d  F,  where  summation 
runs  over  e  hnite  elerhents  of  volume  V_. 

(ii)  Compute  T  c  —  J2e  Id  v  dahKc  with  A  the  external  surface  area  of 

element  e,  with  AT'  the  incremental  traction  vector. 

(iii)  Compute  B  ^  —  jy  AR'g^i,Nl.  d  V,  with  AB“  the  incremental  body  force 

vector.  ~ 

(iv)  Compute  D,.  =  T.e  Iv[^)ha^\i) +JaJA  “  d  Z,  with  y 

_ — 

and  Cf,^  both  found  from  step  1(a). 

(v)  Invert  Eq.  (60)  and  solve  for  Au'^,  i.e. 

Au  =  K-‘(T  +  B  +  D).  (61) 


(c)  Iterate  to  step  1(a)  until  Au„+i  —  Aii„  — ^  0,  where  n  and  n+  1  denote  iterations. 

2.  Update  the  total  nodal  displacement  held  =  u,  +  Au. 

3.  Project  forward  the  intermediate  conhguration  (F**  and  defect  arrangement  |)  for  the  next 
increment  at  each  element  integration  point  (i.e.  unit  cell),  based  on  the  current  local  state. 

4.  Return  to  step  1  until  hnal  load  is  achieved. 

No  problems  are  introduced  in  our  approach  due  to  the  fact  that  the  held  F**  is 
incompatible,  i.e.,  that  in  general.  In  continuum  crystal  plasticity  theories, 

the  mapping  F**  is  generally  not  required  to  be  compatible,  and  in  fact,  the  skew-symmetric 
gradient  of  F**  measuring  its  nonintegrability  is  often  associated  with  geometrically  necessary 
dislocation  density  tensor  (Gurtin,  2002).  In  the  unit  cell  calculations  presented  here,  F**  is 
imposed  uniformly  over  each  cell,  and  F^  includes  the  effects  of  the  hne  scale  perturbation  v, 
except  in  the  hrst  iteration  of  the  numerical  procedure  described  above,  where  the  Born 
assumption  is  used  as  the  initial  guess  for  updating  the  atomic  coordinates.  We  only  consider 
isolated  unit  cells  in  the  applications  that  follow,  but  the  theory  invokes  no  restrictions  on 
how  F**  or  its  gradients  can  vary  from  cell  to  neighboring  cell  on  the  coarse  scale. 


6.  Application:  defects  in  tungsten 

In  the  present  work  our  formulation  is  applied  to  study  the  mechanical  behavior  of  pure 
tungsten  (W),  a  BCC  transition  metal  of  relatively  high  mass  density.  Its  combination  of 
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high  density  and  strength  make  it  an  attractive  material  for  use  in  defense  applications 
such  as  ordnance  (Zhou  and  Clifton,  1997;  Clayton,  2005).  The  potential  energy  function 
describing  atomistic  interactions  in  W  is  discussed  in  Section  6.1.  Numerical  results  for 
illustrative  boundary  value  problems  then  follow.  In  Section  6.2,  these  describe  the 
nonlinear  elastic  response  of  a  deforming  W  crystal  containing  (i)  vacancies,  (ii)  screw 
dislocations  of  like  sign  and  screw  dislocation  dipoles,  and  (hi)  low-angle  twist  boundaries. 
Then  in  Section  6.3  we  address  the  energetic  and  stiffness  properties  of  a  plastically 
deforming  single  crystal  containing  screw  dislocations. 


6.1.  Atomistic  potential:  tungsten 

We  employ  an  empirical  [\l-body  potential  specifically  developed  for  transition  metals 
(Finnis  and  Sinclair,  1984)  in  order  to  estimate  the  free-energy  potential  of  Eq.  (39). 
Duesberry  and  Vitek  (1998)  exercised  the  Finnis-Sinclair  potential  to  model  screw 
dislocation  core  structures,  construct  generalized  stacking  fault  energy  surfaces,  and 
predict  plastic  slip  anisotropy  (i.e.  tension-compression  asymmetry  in  the  yield  function)  in 
multiple  BCC  transition  metals,  including  W.  Liu  et  al.  (2004)  used  the  Finnis-Sinclair 
potential  for  studying  edge  dislocation  glide  in  W  and  molybdenum,  and  Tian  and  Woo 
(2004)  invoked  a  variation  of  this  potential  (Ackland  and  Thetford,  1987)  to  simulate 
screw  dislocation  motion  in  W.  Please  note  that  while  the  choice  of  atomistic  potential  is 
extremely  important  in  order  to  ensure  the  appropriate  physics  are  captured,  our  primary 
interest  in  the  present  research  effort  is  development  of  the  multiscale  homogenization 
technique,  which  may  be  exercised  using  any  available  potential  energy  function  of  form 
(39),  or  more  generally,  for  any  crystalline  material  for  which  atomic  energy,  force,  and 
stiffness  are  known  (e.g.,  in  terms  of  tabulated  values  obtained  from  more  refined 
electronic  structure  calculations).  The  particular  potential  used  here  (Finnis  and  Sinclair, 
1984)  was  selected  because  of  its  apparently  adequate  capability  for  describing  dislocation 
energetics  in  W,  as  reported  in  the  above-mentioned  papers.  The  total  potential  energy  E 
of  a  set  of  atoms  at  positions  for  j  —  1,2, ...,  I\1  is  given  by  the  sum 

E^En  +  Ep  (62) 

where  E^  is  the  I\l-body  term  that  is  a  function  of  the  superposition  of  the  local  electronic 
charge  density  p,  the  latter  obtained  from  a  further  superposition  of  atomic  charge 
densities,  (j).  Also  in  Eq.  (62),  Ep  is  a  pair  potential  that  models  core-core  interactions.  The 
f\l-body  term  is  constructed  as 

En  ^ -Aj2f(P{i)l  (63) 

(j) 


where 

fiPij)) 


P(j)  =  ^(''(/■w) 


(64) 


is  always  non-negative  and  real, 

'"(A*)  ~  \^{i\k)  \  —  |z(/c)  —  Z(/)l, 


(65) 
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0(r)  = 


{r  —  df,  r^d, 
0,  r>d, 


(66) 


and  A  is  an  empirical  constant.  The  adjustable  parameter  d  represents  a  cut-off  zone  for 
superposition  of  local  charge  densities.  For  BCC  tungsten  with  lattice  parameter  a,  the 
value  of  d is  chosen  to  lie  between  the  second-  and  third-nearest  neighbors,  i.e.  a<d< a-Jl, 
for  computational  convenience.  The  pair  potential  Ep  is  formed  as 

(67) 

j,kj*k 

where  is  of  the  polynomial  form 

{(r  — c)^(co-t-Cir-|-C2?'^),  r^c, 

t  ttO  1  2  t, 

0,  r>c, 

with  empirical  constants  Co,  C\,  and  C2-  The  parameter  c,  also  adjustable,  is  presently 
assigned  a  value  between  the  second-  and  third-nearest-neighbor  distances.  Finnis  and 
Sinclair  (1984)  obtained  the  other  constants  via  calibration  to  experimental  data  on 
macroscopic  elastic  properties;  Tables  1  and  2  list  the  experimental  and  fitted  parameters, 
respectively. 

Notice  that  E  is  the  total  energy  of  N  atoms  in  the  fine  scale  unit  cell.  The  Helmholtz  free 
energy  density  in  a  continuum  sense,  T,  is  related  to  E  as 

(69) 

where  U  is  the  continuum  internal  energy,  fj  is  the  continuum  entropy  per  unit  volume,  and 
6  is  the  absolute  thermodynamic  temperature  of  the  system,  which  we  assume  is  zero  in  the 


Table  1 

Experimental  lattice  quantities  for  W 


A 

Lattice  parameter  (A) 

3.1652 

If 

Cohesive  energy  (eV) 

-8.90 

Cu 

Elastic  constant  (GPa) 

522.4 

C\2 

Elastic  constant  (GPa) 

204.4 

C44 

Rhombohedral  shear  modulus  (GPa) 

160.6 

c 

Tetragonal  shear  modulus  (GPa) 

159.0 

B 

Bulk  modulus  (GPa) 

310.4 

Pc 

Cauchy  pressure  =  j(Ci2  —  C44) 

21.9 

Table  2 

Constants  for  EAM  potential  (W) 

d(A) 

4.40024 

A  (eV) 

1.896373 

c(A) 

3.25 

Co 

47.1346499 

C| 

-33.7665655 

C2 

6.2541999 
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last  of  Eq.  (69)  for  the  present  case  of  molecular  statics  simulations,  such  that  NU  —  E. 
Substitution  of  Eq.  (69)  into  earlier  definitions  then  gives 
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to  first  order, 


(72) 


(73) 


Derivatives  of  E  with  respect  to  atomic  displacements  are  listed  by  Finnis  and 
Sinclair  (1984)  and  are  not  repeated  here.  In  our  numerical  implementation,  expressions 
(70)-(73)  are  evaluated  analytically  and  subsequently  used  in  Eqs.  (56)  and  (60). 


6.2.  Multiscale  nonlinear  elasticity  with  defects 

Here  we  report  results  of  simulations  of  unit  cells  containing  various  numbers  of  atoms 
configured  to  represent  several  classes  of  crystal  defects.  Essentially,  we  perform 
calculations  in  which  the  tangent  stiffness  of  a  single  Lagrangian  finite  element  integration 
point  is  determined  by  the  microscopic  (i.e.,  atomistic)  response.  A  sample  atomistic  unit 
cell  representative  of  all  calculations  described  in  the  present  work  is  illustrated  in  Fig.  3. 
The  BCC  unit  cell  is  a  rectangle  of  dimensions  Li  x  L2  x  —  ayf2N\  x  aV6N2  X  aVlNi, 
where  Ni,  N2,  and  A3  are  respectively  the  number  of  repeating  planes  stacked  in  the  [1  1  1]-, 
[1  1  2]-,  and  [1  T  0]-directions.  Periodic  boundary  conditions  are  applied  along  all  faces  of 
the  unit  cell  in  the  usual  manner  (cf.  Vitek,  1976),  such  that  atoms  exiting  the  unit  cell 
during  the  calculation  are  mapped  back  into  the  cell  on  the  opposite  face,  thereby 
preserving  the  total  mass  of  the  system. 

In  the  present  section  (6.2),  we  examine  the  response  of  W  containing  defects  under 
the  conditions  of  null  dislocation  flux,  i.e.  the  local  deformation  gradient  F  =  =  FF  in 

the  context  of  Eq.  (33)  and  Fig.  2.  Our  goal  here  is  to  examine  aspects  of  material 
behavior  that  could  be  used  subsequently  in  stand-alone  continuum  theories,  in  particular 
details  associated  with  stored  energy  of  defect  fields  and  the  effects  of  applied  defor¬ 
mations  on  elasticity  and  stored  energy,  for  various  fixed  defect  concentrations.  Please 
recall  that  our  approach  is  presently  restricted  to  static  calculations,  i.e.,  0  =  0  K, 
as  we  have  yet  to  develop  kinetics  relations  and  scaling  laws  for  finite  temperature 
response.  Nonetheless,  one  may  still  draw  conclusions,  at  least  in  a  qualitative  sense, 
regarding  stress  and  energetics  associated  with  dislocation  behavior  in  BCC  metals 
in  the  context  of  lattice  statics  calculations  (Vitek,  1976;  Duesberry  and  Vitek,  1998).  A 
standard  trend  in  the  literature  has  been  development  of  such  scaling  methods  for  the 
purely  mechanical  problem  before  extending  to  the  finite  temperature  regime  (Dupuy  et  ah. 
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Fig.  3.  Atomistic  scale  unit  cell  for  perfect  BCC  lattice. 


2005).  Milstein  and  co-workers  (Milstein  and  Chantasiriwan,  1998;  Milstein  et  al., 
2004)  have  analyzed  the  response  of  a  number  of  metallic  crystalline  elements  using 
molecular  mechanics  techniques  (including  lattice  statics  with  the  embedded 
atom  method);  nonlinear  elastic  moduli,  void  formation,  and  elastic  instabilities 
were  studied  in  these  works,  though  periodic  initial  defect  distributions  and  dislocations 
were  not. 

In  the  present  set  of  calculations,  initial  atomic  coordinates  are  found  using  a  two-step 
procedure:  first  the  linear-elastic  solution  for  displacement  field  of  the  defect  is  applied  to 
the  atoms,  then  a  conjugate  gradient  algorithm  (Plimpton  and  Hendrickson,  1993)  is  used 
to  transition  the  atomic  positions  to  a  stable  local  minimum  energy  state.  Subsequently, 
the  response  to  applied  deformation  is  computed  using  our  asymptotic  homogenization 
scheme  according  to  the  numerical  procedure  described  in  Section  5.  The  applied  (i.e., 
coarse-scale)  deformation  gradient  field  (in  conjunction  with  fine  scale  periodicity)  is 
uniaxial  stretching  in  the  [1  1  l]-direction  (see  Fig.  3)  over  a  range  of  l.OOO^Tii  ^  1.025, 
with  the  lateral  edges  fixed  (covariant  Cartesian  notation  is  used  here  and  in  subsequent 
figures  for  simplicity,  i.e.  F  y  — ^  Tn).  We  also  compare  the  final  configurations  attained 
using  our  procedure  with  that  of  incremental  energy  minimization.  In  the  latter  approach, 
a  small  increment  in  the  stretch  field  is  first  imposed  uniformly  over  all  atoms,  and  then  a 
conjugate  gradient  program  is  used  to  update  the  atomic  coordinates  to  the  corresponding 
local  minimum  energy  state.  This  process  continues  with  a  new  set  of  conjugate  gradient 
minimization  iterations  conducted  upon  application  of  each  successive  stretch  increment 
until  the  final,  fully  deformed  configuration  is  attained. 

First  we  consider  a  point  vacancy,  as  shown  in  Fig.  4(a).  The  initial  configuration  is 
constructed  simply  by  removing  the  atom  closest  to  the  centroid  of  the  unit  cell.  The  defect 
density  in  this  case  is  defined  as  the  volume  fraction  of  missing  atoms,  i.e.,  p'*  =  1/f^, 
where  I\l  is  the  total  number  of  atoms  prior  to  vacancy  creation,  ranging  from  2016  to 
32256  among  the  simulations  discussed  here.  Defect  energy  is  shown  in  Fig.  4(b),  defined 
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Fig.  4.  Vacancy  configuration  (a),  strain  energy  comparison  between  computational  methods  (b),  defect  energy 
(c),  elastic  stiffness  (d),  and  Zener  anisotropy  (e). 


on  a  per-atom  basis  as 


E-E 


(74) 


where  E  is  the  total  potential  energy  of  the  system  and  E  is  the  total  potential  energy  of  a 
perfect  BCC  W  lattice  of  the  same  dimensions  and  same  number  of  atoms  (prior  to 
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vacancy  formation),  subjected  to  the  same  deformation  boundary  conditions.  From 
Fig.  4(b),  we  see  that  our  homogenization  approach  is  validated  in  the  sense  that  it  predicts 
minimum  energy  atomic  configurations  that  compare  favorably  with  those  obtained  using 
incremental  conjugate  gradient  minimization.  In  Fig.  4(c),  the  total  defect  energy  is 
shown.  Schultz  (1991)  reported  an  experimental  value  of  3.6  eV  for  vacancy  formation 
energy  in  W  at  null  applied  strain.  Notice  that  the  defect  energy  increases  ^yith  applied 
stretch  T^ii.  Shown  in  Fig.  4(d)  is  the  elastic  modulus  component  C""  =  Cjj,  computed 
via  Eq.  (72).  Values  decrease  drastically  with  applied  stretch,  to  around  80%  of  their 
original  magnitudes  at  Fn  =  1.025.  In  contrast,  the  defect  density  has  little  effect  on  C*"^ 
though  this  could  be  expected  at  such  dilute  vacancy  concentrations  as  are  considered  here 
(i.e.,  porosities  less  than  0.1%).  Fig.  4(c)  shows  the  Zener  anisotropy  factor  A,  dehned  by 
(cf.  Hirth  and  Lothe,  1982) 


A  = 


2C 


1212 


2C, 


-22 


-  c 


1122 


Cn 


a12  ■ 

Ci2 


(75) 


where  the  Cartesian  contravariant  notation  is  used  in  Fig.  4(c)  and  in  the  expression 
following  the  first  equality  in  Eq.  (75).  Note  that  defect-free  W  is  nominally  isotropic  at 
null  deformation,  i.e.,  A  —  1.00,  and  that  the  anisotropy  increases  dramatically  as  the 
lattice  is  stretched.  This  result  may  in  part  be  due  to  a  limitation  in  the  Finnis  and  Sinclair 
(1984)  potential,  whose  parameters  are  calibrated  to  small-strain  elastic  moduli,  although 
first-principles  (i.e.,  quantum-mechanical)  calculations  have  indicated  a  departure  from 
isotropy  in  W  at  large  pressures  (Ruoff  et  ah,  1998). 

Regarding  the  vacancy  concentrations  considered  here  (p‘*  on  the  order  of  lO”"*),  we  note 
that  such  high  concentrations  would  be  unlikely  to  exist  in  commercially  melt-grown  pure 
single  crystals,  and  that  for  very  dilute  concentrations  (on  the  order  of  parts  per  million), 
the  effects  on  the  energetic  and  mechanical  properties  would  apparently  be  quite  small. 
However,  for  tungsten  crystals  and  alloys  generated  from  powder  consolidation/ 
compaction  and  sintering  processes,  even  higher  porosities  on  the  order  of  several  percent 
are  not  uncommon  (cf.  Yih  and  Wang,  1979).  In  such  cases,  the  effects  on  properties  would 
be  substantial,  perhaps  orders  of  magnitude  more  so  than  the  results  listed  here  in  Fig.  4. 
Such  highly  porous  tungsten  is  used,  for  example,  in  electrical  applications  (Selcuk  and 
Wood,  2005),  where  cyclic  fatigue  is  of  concern. 

The  response  of  W  containing  a  periodic  array  of  screw  dislocations  is  investigated  next. 
The  defect  conhguration  is  illustrated  in  Fig.  5(a).  The  dislocation  tangent  line  and  burgers 
vector  b  are  oriented  along  the  [1  1  l]-direction  and  pass  through  the  centroid  of  the  unit 
cell,  the  latter  having  a  magnitude  of  h  =  |b|  =  \/3a/2.  Initial  atomic  positions  are 
prescribed  via  the  usual  displacement  held  solution  for  a  screw  dislocation  embedded  in  an 
inhnite  isotropic  elastic  body,  i.e.  u  —  bO/ln,  where  6  is  an  angular  coordinate  about  the 
axis  of  the  dislocation  line.  The  scalar  dislocation  density  is  dehned  as  the  defect  line  length 
per  unit  reference  volume,  i.e.,  p'*  =  The  defect  energy  per  atom,  dehned  as  in 

Eq.  (74)  and  shown  in  Fig.  5(b),  is  computed  accurately  by  our  homogenization  scheme  as 
is  verihed  by  the  incremental  conjugate  gradient  solutions.  Furthermore,  we  see  a  linear 
increase  in  stored  defect  energy  with  applied  deformation  Fn,  and  a  roughly  linear  increase 
in  F^  with  increasing  dislocation  density  p‘*.  Gibeling  and  Nix  (1980)  discussed,  from  the 
standpoint  of  discrete  dislocation  modeling,  how  the  strain  energy  supported  by 
dislocations  may  be  amplihed  by  applied  external  deformations,  and  Clayton  (2005) 
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Fig.  5.  Screw  dislocation  configuration  (a),  strain  energy  comparison  between  computational  methods  (b),  defect 
energy  (c),  elastic  stiffness  (d),  and  Zener  anisotropy  (e). 


assumed  a  linear  dependence  of  stored  elastic  energy  on  dislocation  density  in  a  continuum 
crystal  plasticity  model  of  single  crystalline  W.  Note  that  we  are  not  attempting  to  compute 
isolated  dislocation  energies  as  has  been  the  goal  of  previous  studies  (cf.  Cai  et  ah,  2003). 
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Instead,  here  we  wish  to  include  the  effects  of  neighboring  defects,  thereby  representing  a 
distribution  of  lattice  imperfections  as  would  occur  in  a  plastically  deforming  sample  of 
material.  Note  also  that  a  distribution  of  identically  oriented  dislocations  produces  a  net 
(geometrically  necessary)  dislocation  density  tensor  in  the  sense  of  Nye  (1953),  implying 
curvature  in  the  lattice  (see  also  Ashby,  1970).  From  Figs.  5(d)  and  (e),  we  see  that  the 
dislocation  density  tends  to  accelerate  the  decrement  in  stiffness  value  C'  * "  with  increasing 
applied  stretch,  whereas  the  anisotropy  A  tends  to  be  suppressed  by  increasing  dislocation 
content  as  the  stretch  increases.  In  particular,  at  null  applied  strain,  C'"*  decreases  from 
514  to  500  GPa  as  the  defect  density  is  increased  from  0  to  0.026/nm^.  At  an  applied  stretch 
of  Fii  =  1.025,  C*"'  decreases  from  414  to  395  GPa  over  this  same  range  of  Although 
the  trend  of  decreasing  elastic  modulus  with  dislocation  content  has  been  reported 
elsewhere  in  the  literature  following  physical  experiments  (Smith,  1953)  and  analytical 
continuum  modeling  (Lebedev,  1996),  it  has  not  been  emphasized  previously  in  continuum 
plasticity  models.  It  is  also  noted  that  the  maximum  elastic  deformation  attained  in  these 
simulations,  i.e.,  2.5%  stretch,  is  larger  than  would  be  expected  in  real  crystals  containing 
mobile  defects  wherein  yielding  and  plastic  flow  would  occur  more  readily.  Atomic 
motions  are  constrained  here  due  to  the  null  temperature  prescription  and  uniaxial  strain 
boundary  conditions,  though  our  results  are  consistent  with  those  of  other  researchers 
under  similar  constraints  (cf.  Vitek,  1976). 

Results  for  a  unit  cell  containing  a  screw  dislocation  dipole  are  presented  now.  The 
defect  configuration  is  shown  in  Fig.  6(a).  The  dislocation  tangent  lines  are  oriented  along 
the  [1  1  l]-direction.  The  first  dislocation  of  the  pair  is  located  at  (I/4L2, 1/4L3),  and  the 
second  is  located  at  (3/41,2,3/41,3),  thus  maintaining  a  minimum  separation  distance  of 
(I/2L2, 1/21,3)  between  the  dislocations  located  within  the  cell  as  well  as  between  image 
dislocations  implied  by  the  periodic  boundaries.  The  atomic  displacements  are  initially 
prescribed  as  a  superposition  of  the  isotropic  linear-elastic  displacement  held  solutions  of 
the  two  defects,  with  the  burgers  vector  of  the  first  dislocation  oriented  positively  in  the 
[1  1  l]-direction  and  that  of  the  second  oriented  negatively  in  the  [1  1  l]-direction.  In  this 
case,  the  dislocation  density  p'*  =  2/(Z,2Z,3),  and  the  net  Nye  (1953)  tensor  vanishes,  as  the 
burgers  vectors  cancel  out,  meaning  the  dislocation  density  is  “statistically  stored”  in  the 
sense  of  Ashby  (1970).  The  defect  energy  computed  with  our  homogenization  scheme  is 
plotted  in  Fig.  6(b);  again,  the  numerical  results  are  validated  by  conjugate  gradient 
molecular  statics  simulations.  From  Fig.  6(c),  we  see  that  the  defect  energy  per  atom 
increases  roughly  linearly  with  the  dislocation  density,  p‘*,  and  it  increases  linearly  with 
applied  stretch  Fi  i .  The  ratio  of  defect  energy  to  dislocation  density  is  naturally  smaller  for 
the  dislocation  dipole  than  for  the  single  screw  dislocation  (Fig.  5),  as  the  local  stress  fields 
of  the  dislocations  comprising  the  dipole  counteract  to  a  certain  degree.  For  example,  at 
Fii  —  1.025,  for  a  single  dislocation,  F^  —  0.0254 eV/atom  at  p'*  =  0.026/nm^,  whereas  for 
a  dipole,  the  energy  is  only  77*  =  0.0214  eV/atom  at  twice  the  density  p'*  =  0.052/nm^.  As 
shown  in  Fig.  6(d),  the  stiffness  coefficient  C"'*  decreases  with  increasing  defect  density 
and  increasing  stretch.  Anisotropy  A  (Fig.  6(d))  increases  with  stretch  and  increases  with 
dislocation  content  at  low  values  of  p'*  (^0.013/nm^),  but  decreases  as  p'*  is  increased 
further. 

The  final  defect  configuration  examined  here  is  a  low  angle  twist  grain  boundary.  In 
Fig.  7(a),  we  describe  this  boundary  type  using  the  disclination  concept  (Li,  1972;  Clayton 
et  ah,  2006).  Considered  here  are  twist  disclinations,  the  rotational  analogs  of  screw 
dislocations.  The  defect  engenders  a  misorientation  across  the  (1  1  0)-plane  of  strength 
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Fig.  6.  Screw  dislocation  dipole  configuration  (a),  strain  energy  comparison  between  computational  methods  (b), 
defect  energy  (c),  elastic  stiffness  (d),  and  Zener  anisotropy  (e). 


CO  =  0.2  rad,  the  value  of  which  was  chosen  for  convenience  such  that  a  stable  local 
minimum  energy  configuration  could  be  found  for  the  initial  atomic  arrangement.  In  the 
disclination  description,  the  Frank  vector  is  oriented  parallel  to  the  [1  1  0]-direction,  the 
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Fig.  7.  Twist  disclination  configuration  (a),  strain  energy  comparison  between  computational  methods  (b),  defect 
energy  (c),  elastic  stiffness  (d),  and  Zener  anisotropy  (e). 


disclination  line  i  is  oriented  along  the  [1  1  l]-direction,  and  the  axis  of  twist  passes  through 
the  centroid  of  the  unit  cell.  The  disclination  line  density  is  simply  =  \/{L2L^).  In  our 
simulations,  the  atoms  are  initially  displaced  according  to  the  elastic  solution  (Fig.  7(a)), 
then  subjected  to  an  initial  energy  minimization  via  a  conjugate  gradient  solver.  As  the 
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deformation  proceeds,  the  atomic  coordinates  are  updated  via  the  homogenization  scheme 
whose  algorithm  is  outlined  in  Section  5.  The  defect  energy  results  are  again  validated 
via  comparison  with  an  incremental  conjugate  gradient  solution  in  Fig.  7(b).  We  see  in 
Fig.  7(c)  that  the  defect  energy  per  atom  increases  nonlinearly  with  disclination  line 
density;  in  other  words,  it  increases  nonlinearly  with  grain  boundary  area  per  unit  volume. 
Note  that  we  are  essentially  simulating  a  periodic  array  of  such  low-angle  grain 
boundaries,  as  have  been  known  to  appear  in  pure  W  subjected  to  severe  plastic 
deformation  (Valiev  et  ah,  2002;  Wei  et  ah,  2005,  2006).  In  contrast  to  the  other  classes  of 
defects  examined  in  the  present  work,  the  defect  energy  decreases  with  applied  stretch,  and 
the  stiffness  component  increases  with  increasing  defect  density.  From  Fig.  7(d),  at 
null  applied  strain,  C" '  ’  increases  from  5 14  to  520  GPa  as  the  defect  density  is  increased 
from  0  to  0.026/nm^.  At  an  applied  stretch  of  Fn  =  1.025,  C""  increases  over  the  same 
range  of  p‘^  from  414  to  419  GPa.  Anisotropy,  shown  in  Fig.  7(e),  is  suppressed  slightly 
with  increasing  disclination  content. 

6.3.  Multiscale  elastoplasticity 

Next  we  model  the  response  of  screw  dislocations  in  the  context  of  elastoplasticity.  The 
full  decomposition  F  =  =  FF  of  Eq.  (33)  and  Fig.  2  applies.  Initial  conditions  for  the 

unit  cell  are  identical  to  those  of  Fig.  5(a):  a  single  screw  dislocation  is  oriented  in  the 
[1  1  l]-direction,  the  atoms  are  initially  displaced  according  to  the  isotropic  elastic  solution 
for  the  displacement  field  of  a  screw  dislocation,  and  then  subjected  to  conjugate  gradient 
minimization  at  null  applied  macroscopic  strain.  Again,  the  dislocation  line  density  is 
found  as  =  \/{L2L^).  Here  we  apply  pure  shear  deformation  in  the  1-3  plane,  over  the 
range  0.00  ^0.10.  Periodic  boundary  conditions  are  again  invoked  across  all  faces  of 

the  unit  cell  in  the  atomistic  domain.  For  the  case  of  pure  shear  in  the  1-3  plane,  the 
multiplicative  decomposition  F  =  F^F**  yields 

-Fi3  =  y  =  +  y^,  (76) 

where  is  the  lattice  shear  associated  with  the  external  stress  and  y^  is  the  cumulative 
plastic  slip,  assumed  in  our  idealized  problem  here  to  occur  as  a  result  of  dislocation  glide 
on  (1  1  0)-planes.  Plastic  slip  is  assumed  to  follow  the  simple  kinetic  relation 

,  f  0  (Vy^y^), 

y  =  <  .3  .  (77) 

I  7  -  7o  -  ^  1  +  sin[(27i/ ab)(y  -  y^  -  a6/4)]}  (Vy  ^ y^), 

where  y^  is  the  initial  yield  strain,  yj  is  an  additional  lattice  strain  required  to  overcome  the 
Peierls  barrier,  and  a  =  1  /F3  scales  the  amount  of  shear  strain  accumulated  in  the  unit  cell 
due  to  the  slip  of  a  dislocation  over  the  distance  of  one  burgers  vector.  As  will  be  shown 
later,  Eq.  (77)  results  in  a  shear  stress-shear  strain  relationship  that  oscillates  in  a 
sinusoidal  manner,  in  accordance  with  the  Peierls  model  of  periodic  lattice  resistance  stress 
discussed  by  Hirth  and  Lothe  (1982).  No  rate  dependence  is  included  in  Eq.  (77),  which  is 
restricted  to  monotonic  shear.  We  regard  Eq.  (77)  as  a  minimal  substitute  for  dislocation 
kinetics,  which  are  not  accounted  for  explicitly  in  our  present  multiscale  scheme,  the  latter 
limited  to  0  K  calculations  at  the  atomic  scale.  However,  such  an  approximation  is  deemed 
adequate  to  study  strain  energy  and  associated  defect  properties  in  presence  of  plastic 
deformation,  and  not  plasticity  kinetics  itself  which  may  require  molecular  dynamics  at  the 
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fine  scale  in  our  framework.  The  simulations  discussed  here  assume  that  the  primary 
dislocation  remains  hxed  in  the  unit  cell  (in  the  atomistic  domain),  and  that  no  defect 
generation  or  annihilation  occurs,  as  other  screw  dislocations  convect  through  the  lattice 
and  cause  an  increase  in  plastic  shear  In  other  words,  time  differentiation  of  Eq.  (77) 
produces  the  evolution  equation  for  the  plastic  deformation  gradient  F*”  (54)  and 
corresponding  intermediate  conhguration  update  (Fig.  2),  and  —  0  may  be  regarded  as 
part  of  the  general  set  of  evolution  equations  (55).  Here,  following  the  molecular  statics 
calculations  of  Vitek  (1976),  we  choose  —  0.020  and  yj  =  0.009,  the  former  a  shear  strain 
required  to  cause  [1  1  1](1  1  0)-screw  dislocation  core  shuffling  in  BCC  W,  the  latter 
following  his  result  that  a  total  shear  strain  of  0.029  is  required  to  sustain  steady 
dislocation  glide. 

The  1-3  shear  component  of  the  first  Piola-Kirchhoff  stress,  is  shown  in  Fig.  8(a), 
where  components  of  the  stress  tensor  P  are  computed  from  Eq.  (24)  and  (69)  as 


Pd  =  J 


P-\  rP-lA 


0F; 


La  .0! 


jP-\  1  +  Ep)  f_i^ 


(78) 


Prior  to  initial  yield,  slight  differences  in  slope  of  the  stress-strain  response  are  evident 
among  the  three  curves,  with  each  curve  corresponding  to  a  different  dislocation  density  p^. 


(a)  F 13  (b)  F (3 


W  F 13 


Fig.  8.  Results  from  elastic-plastic  homogenization  simulations:  shear  stress  (a),  elastic  shear  stiffness  (b),  and 
total  relative  elastic  energy  (c). 
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Such  differences  are  due  in  part  to  the  decrease  in  effective  elastic  shear  modulus,  = 

C33 ,  with  increasing  defect  density,  as  is  seen  in  Fig.  8(b).  Upon  initial  yield,  the  stress-strain 
curves  follow  Peierls-type  oscillations  (a  direct  result  of  imposition  of  Eq.  (77)),  with  each 
period  corresponding  to  the  passage  of  a  single  [1  1  l]-screw  dislocation  across  the  unit  cell, 
gliding  on  a  (1  1  0)-plane.  As  the  dislocation  density  p‘^  is  inversely  related  to  the  dimension 
L3  of  the  unit  cell  in  our  calculations,  the  larger  the  dislocation  density,  the  smaller  the  unit 
cell  size,  the  larger  increment  in  shear  strain  associated  with  the  passage  of  each 
dislocation,  and  the  longer  the  period  of  shear  stress  oscillations  in  Fig.  8(a).  In  the  absence 
of  such  oscillations,  the  material  behavior  would  be  virtually  elastic-perfectly  plastic.  The 
total  energy  per  atom,  —  E/N  —  U'^,  is  shown  in  Fig.  8(c),  where  the  cohesive  energy 
U‘^  —  — 8.90eV/atom  (Table  1)  has  been  subtracted  such  that  a  defect-free,  undeformed 
lattice  would  have  null  energy  F*.  At  zero  applied  shear,  differences  among  energy  curves  are 
due  purely  to  differences  in  defect  concentration,  while  the  energy  oscillations  at  larger  shear 
deformations,  Fi3>  0.020,  are  associated  with  the  lattice  deformation  (F^)  appearing  in 
conjunction  with  the  applied  stress. 

7.  Discussion 

We  discuss  here  the  relevance  and  potential  utility  of  our  numerical  results,  specifically 
the  properties  computed  and  presented  above  for  tungsten  crystals  with  defects.  We  then 
review  the  advantages  and  differences  of  our  method  in  comparison  to  existing  multiscale 
approaches. 

The  fundamental  results  we  have  presented  here  describe  the  thermodynamic  free  energy 
and  its  second  derivative  (i.e.,  elastic  stiffness)  of  unit  cells  of  BCC  tungsten  containing 
various  point,  line,  and  area  defects.  In  contrast,  in  the  continuum  plasticity  literature,  the 
primary  area  of  research  has  historically  been  the  plastic  flow  rule  and  associated 
quantities  (e.g.,  yield  surface,  strain  hardening  parameters,  and  internal  variables  linked  to 
these),  as  the  flow  rule  dominates  the  stress-strain  response  when  the  elastic  strains  are 
small,  as  is  usually  the  case  for  engineering  metals  undergoing  even  finite  quasi-static 
deformations.  Our  method  has  not  been  applied  towards  determination  of  plastic  kinetic 
relationships  and  parameters  such  as  these.  Nonetheless,  the  residual  energy  of  defects  that 
we  study  is  directly  related  to  the  stored  energy  of  cold  working  and  the  fraction  of  stress 
power  converted  to  heat  energy,  an  important  aspect  of  continuum  plasticity  modeling  (see 
e.g.  Rosakis  et  ah,  2000).  Fundamental  relationships  between  stored  energy  of  cold  work 
and  microstructure  properties  such  as  dislocation  content  are  not  readily  available  in  the 
existing  literature  for  most  materials,  and  our  method  and  results  offer  a  way  to  deduce 
such  relationships.  Furthermore,  in  strain  gradient-type  crystal  plasticity  models  (see  e.g., 
Gurtin,  2002;  Clayton  et  al.,  2006),  derivatives  of  the  free  energy  with  respect  to  the 
(geometrically  necessary)  dislocation  density  (i.e.,  thermodynamic  conjugate  forces)  are 
assumed  to  give  rise  to  internal  stresses  and  in  some  instances,  directional  strain  hardening 
with  a  size  or  length  scale  dependence.  Our  technique  is  conceptually  able  to  compute  the 
energy  that  would  give  rise  to  such  an  effect,  though  we  do  not  explicitly  compute  such 
conjugate  forces  here  in  the  context  of  idealized,  perfectly  periodic  defect  configurations. 

The  dislocation  densities  simulated  here  are  large  (e.g.,  on  the  order  of  10*®/m^)  in 
comparison  to  contents  found  in  homogeneously  deforming  single  crystals.  For  example. 
Argon  and  Maloof  (1966)  report  values  of  on  the  order  of  10*^/m^  for  pure  tungsten 
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single  crystals  deformed  up  to  5%  axial  tensile  strain.  Our  simulated  values  correspond  to 
dislocation  spacings  on  the  order  of  lOnm.  However,  we  argue  that  in  highly  strained 
regions  of  crystal,  such  as  in  the  vicinity  of  grain  or  subgrain  boundaries  formed  during 
severe  plastic  deformation  (SPD),  such  dislocation  spacings  are  not  unreasonable.  For 
example,  if  one  considers  a  boundary  comprised  of  sequence  of  dislocations  of  spacing  h  and 
Burgers  vector  b,  the  misorientation  at  the  boundary  can  be  computed  as  bjh  (cf.  Hughes 
et  al.,  2003),  on  the  order  of  4°  for  lOnm-spaced  dislocations  in  tungsten.  Misorientations  of 
such  magnitude  have  been  observed  and  documented  for  subgrain  boundaries  produced  in 
tungsten  crystals  deformed  through  SPD  processes  (Valiev  et  ah,  2002;  Wei  et  ah,  2006). 
Furthermore,  strain  gradient  theories  as  mentioned  above  are  designed  to  address  such 
phenomena,  so  results  presented  here  may  be  applied  to  motivate  continuum  energy 
dependencies,  at  least  in  a  qualitative  sense,  upon  defect  densities  that  serve  an  important 
role  in  such  theories.  Stored  energy  of  cold  working  relative  to  heat  dissipation  may  also 
influence  shear  localization  processes  in  ultra-fine  grained  tungsten  (Wei  et  ah,  2005,  2006), 
which  can  exhibit  dislocation  densities  of  the  magnitude  studied  here. 

Our  theory  has  no  limitations  with  regards  to  defect  density,  and  conceptually  permits 
the  consideration  of  inhnitely  lower  defect  densities  by  increasing  the  number  of  atoms  in 
the  unit  cell.  However,  computational  expense  associated  with  the  present  numerical 
implementation  prevents  us  from  considering  problem  sizes  on  the  order  of  tens  of  millions 
of  atoms,  as  would  be  needed  to  reach  very  low  defect  densities  typically  observed  in 
commercially  pure  melt-grown  single  crystalline  tungsten  (Argon  and  Maloof,  1966).  It 
should  be  noted  that  size  limitations  do  exist  for  alternative  techniques  such  as  molecular 
mechanics,  though  these  more  mature  methods  have  benehted  from  decades  of  research  on 
(parallel)  algorithm  development  and  hence  enjoy  perhaps  greater  computational 
efficiency,  whereas  our  theory,  being  relatively  new,  has  not  yet  been  fully  optimized 
from  a  computational  standpoint. 

Plotted  in  the  results  Sections  6.2  and  6.3  (Figs.  5-8)  are  defect  energies  per  unit 
reference  volume.  This  is  equivalent  to  presentation  of  the  energy  per  atom  if  each  atom 
occupies  a  fixed  reference  volume.  According  to  linear  elasticity  theory,  for  one  screw 
dislocation,  there  is  a  logarithmic  decay  with  distance  from  the  dislocation  line  of  energy 
(per  unit  dislocation  line  length)  due  to  the  stress  field  of  the  defect.  Furthermore,  the 
dislocation  density  in  our  simulations  is  controlled  by  the  dimensions  of  the  unit  cell  (i.e., 
the  bounding  box  of  the  atoms).  For  an  infinitely  periodic  array  of  straight  screw 
dislocations,  the  energy  per  atom  should  not  depend  on  the  system  size  for  a  fixed 
dislocation  density,  so  long  as  the  orientation  of  the  bounding  box  with  respect  to  the 
dislocation  line(s)  and  the  periodicity  and  spacing  of  the  dislocations  are  maintained.  For 
example,  quadrupling  the  area  of  the  box  face  normal  to  the  dislocation  line  and  adding 
three  more  equivalent  dislocations  (for  a  total  of  four  dislocations)  would  give  the  same 
energy  per  atom  as  the  original  box  with  one  dislocation  so  long  as  the  atoms  surrounding 
each  of  the  four  dislocations  in  the  enlarged  box  were  to  occupy  a  perfectly  periodic 
arrangement.  The  energy  would  of  course  depend  on  the  relative  dimensions  (i.e.,  aspect 
ratio)  of  the  edges  of  the  bounding  box;  as  such,  aspect  ratios  were  kept  fairly  constant  in 
the  simulations  in  which  the  absolute  size  of  the  unit  cell  was  used  to  control  the  defect 
density.  Thus,  it  is  noted  that  the  shape  of  the  unit  cell  will  have  an  effect  on  the  defect 
density-versus-strain  energy  relationships,  and  furthermore,  that  size  effect  issues  would 
appear  more  prominently  for  more  complicated  defect  arrangements,  such  as  edge  or 
mixed  dislocations.  We  also  note  that  even  if  such  size  effects  are  unavoidable,  they  are  not 
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necessarily  unnatural,  and  the  choice  of  unit  cell  dimensions  (i.e.,  number  of  atoms)  is  not 
restricted  by  our  method  in  theory,  but  in  practice  due  to  computational  convenience. 

The  effects  of  defects  on  elastic  modulus  and  energy  reported  here,  while  not  profound, 
are  noticeable.  For  example,  the  highest  density  of  dislocations  engenders  a  change  in 
elastic  stiffness  constant  C*”*  on  the  order  of  3%  for  the  unit  cells  deformed  in  tension 
(Fig.  5(d))  and  a  change  in  the  anisotropy  factor  on  the  order  of  10%  (Fig.  5(e)).  For  the 
elastic-plastic  shear  simulations,  differences  in  on  the  order  of  15%  emerge  among 
the  cases  considered  in  Fig.  8(b),  though  these  differences  are  not  highly  evident  in  the 
stress-strain  curves  of  Fig.  8(a),  whose  forms  are  controlled  primarily  through  the  plastic 
flow  prescription  (77).  Flowever,  in  a  load-unload  scenario  (e.g.,  cyclic  fatigue),  the 
changes  in  elastic  moduli  would  conceivably  be  more  apparent  in  the  (unloading  portion  of 
the)  stress-strain  curves.  Furthermore,  the  trends  in  results  presented  here  do  agree 
qualitatively  with  the  experimental  observation  of  decreasing  modulus  with  increasing  cold 
work  and  dislocation  content  reported  elsewhere,  for  other  metals  (Smith,  1953;  Lebedev, 
1996).  The  defect  energies,  while  perhaps  relatively  small  in  comparison  with  the  total 
energy  of  the  system,  are  thought  to  be  of  importance  with  regards  to  continuum  modeling 
of  the  energy  of  cold  working  and  internal  stresses  and  size-dependent  hardening  in  strain 
gradient  crystal  plasticity  models,  as  discussed  above. 

The  initial  state  of  the  material  in  our  simulations  is  one  containing  defects,  though  these 
comprise  a  local  minimum  energy  conhguration  in  the  atomistic  sense.  Relative  to  a  perfect 
defect-free  lattice,  the  system  with  the  defect  contains  positive  total  energy,  though  the 
strain  energies  of  both  the  defective  and  defect-free  systems  are  considered  null  in  the 
undeformed  state.  In  principle,  our  approach  differs  from  many  continuum  plasticity 
implementations,  wherein  the  defect  density  is  presumed  zero  initially  then  and  evolves 
with  F*”,  for  example,  as  the  material  work  hardens.  Flowever,  we  argue  that  since  all  real 
materials  contain  defects  initially,  our  approach  is  realistic  in  this  regard.  The  issue  could 
be  resolved  in  the  context  of  existing  continuum  theories  too,  if  one  assumes  that  the  initial 
yield  stress  is  related  to  the  initial  defect  density,  and  that  the  reference  state  is  offset  by 
some  energy  value  associated  with  this  defect  content.  In  the  context  of  hyperelasticity, 
only  the  derivatives  of  the  energy  with  respect  to  elastic  deformation  contribute  to  the 
stress  state,  and  the  absolute  value  of  the  energy  will  not  have  an  effect.  Though  the  choice 
of  conhguration  may  affect  the  values  of  the  nonlinear-elastic  stiffness  tensor,  it  is 
customary  in  hnite  elastic-plastic  theories  to  express  the  elastic  moduli  in  the  intermediate 
frame  (cf.  Gurtin,  2002).  The  numerical  results  we  present  for  variable  elastic  stiffness 
under  applied  hnite  deformations  are  consistent  with  our  theoretical  formulation,  though 
comparing  our  results  with  published  values  may  not  be  straightforward.  For  example, 
nonlinear  elastic  effects  are  often  computed  experimentally  from  stress  changes  upon 
perturbations  of  the  displacement  held  in  the  spatial  frame  and  then  expressed  in  terms  of 
derivatives  (with  respect  to  strain)  of  the  small  strain  linear  elastic  constants  (Toupin  and 
Rivlin,  1960). 

This  research  effort  has  provided  a  method  for  incorporation  of  the  (continuum) 
kinematics  of  hnite  plasticity  in  conjunction  with  defect  energy  and  nonlinear  elastic 
properties  obtained  directly  from  asymptotic  homogenization  of  the  atomic  scale  response. 
We  do  not  attempt  here  to  derive  continuum  plasticity  kinetics  (e.g.,  how  rules,  yield 
surfaces,  dislocation  velocities  and  generation  rates)  from  the  atomic  scale.  Our  approach 
is  presently  limited  to  static  problems  at  both  hne  and  coarse  length  scales.  The  starting 
point  for  extending  our  method  to  dynamic  problems  would  be  inclusion  of  acceleration  in 
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the  momentum  balance  (Eq.  (4)),  and  following  through  with  the  (nontrivial)  derivation  of 
the  asymptotic  equations  including  a  term  accounting  for  hne  scale  (i.e.,  atomic)  inertia. 
Conceptually,  this  approach  may  enable  one  to  consider  temperature  effects  that  would  be 
essential  for  describing  thermally  activated  dislocation  glide  and  climb.  Additionally,  the 
unit  cell  calculations  would  need  to  be  extended  to  more  complicated  scenarios  (e.g.,  Frank 
Read  sources,  dislocation  pile-ups)  to  capture  strain  hardening  and  dislocation  interactions 
thought  to  be  of  high  importance  in  continuum  plasticity.  The  analysis  discussed  in  Section 
6.3,  dealing  with  the  oscillatory  stress-strain  response,  demonstrates  of  how  our  method 
may  be  applied  in  a  limited  sense  to  plastic  phenomena  without  addressing  the  above 
issues.  A  simple  sinusoidal  plastic  strain  was  prescribed  in  this  case  (77)  as  it  is  perhaps  the 
most  basic,  yet  still  physically  realistic,  flow  model  enabling  us  to  demonstrate  the  energy 
and  stiffness  changes  with  defects  (in  the  context  of  our  multiscale  scheme)  simultaneously 
incorporating  plasticity  kinematics.  The  Peierls-like  oscillatory  stress-strain  curves 
(Fig.  8(a))  are  a  direct  outcome  of  this  choice  of  flow  rule,  though  the  key  results  for 
stiffness  and  energy  (Figs.  8(b)  and  (c))  arise  as  a  consequence  of  the  fine  scale,  atomic 
response  and  would  not  be  available  from  a  purely  continuum-scale  analysis.  A  more 
complex  flow  rule,  perhaps  accounting  for  dislocation  interactions  and  temperature,  could 
easily  be  incorporated  (see  e.g.  Groma  and  Pawley,  1993,  in  the  context  of  discrete 
dislocation  plasticity),  though  it  is  not  clear  that  such  an  exercise  would  provide  additional 
physical  insight  in  the  context  of  the  present  set  of  unit  cell  calculations.  Finally,  we  note 
that  other  multiscale  methods  such  as  the  quasi-continuum  (Shenoy  et  ah,  1999)  or 
bridging  scales  method  (Park  et  al.,  2005)  have  not  been  applied,  at  least  to  our  knowledge 
in  a  universally  accepted  manner,  to  incorporate  finite  plasticity  kinematics  and  kinetics, 
though  it  may  be  quite  possible  for  one  to  do  so. 

Next,  we  further  consider  the  advantages  of  our  method  in  comparison  with  other 
multiscale  modeling  techniques.  The  asymptotic  homogenization  formalism  is  a 
mathematical  method  for  embedding  smaller  scale  information  into  a  larger  scale  (we 
use  “embedding”  here  in  a  materials  modeling  sense,  not  a  mathematical  sense). 
Convergence  properties  related  to  two-scale  continuum  homogenization  have  been 
reported  extensively;  for  example,  Bensoussan  et  al.  (1978)  discusses  how,  in  the  limit 
that  the  scaling  parameter  is  driven  to  zero,  the  asymptotic  solution  converges  to  the  exact 
continuum  solution  for  the  response  (e.g.  local  displacement  field)  of  an  infinitely  periodic 
medium  weakly.  We  also  find  it  important  to  point  out  a  feature  that  distinguishes  our 
work  from  other  atomistic-continuum  multiscale  methods.  That  is  that  our  method  does 
not  require  the  careful  matching  of  displacement  boundary  conditions  between  atomistic 
and  continuum  domains,  and,  instead,  maps  derivatives  (i.e.,  deformation  gradients) 
across  length  scales.  In  our  technique,  the  atomic  coordinates  are  not  required  to  share  the 
same  coordinate  system  and  characteristic  length  scale  as  the  larger  continuum  since 
displacement  is  not  the  key  inter-scale  variable.  Instead,  a  scaling  factor  (2)  here  relates  fine 
and  coarse  scales.  The  coarse  scale  deformation  gradient,  F  is  transmitted  to  the  fine  scale, 
and  extra  degrees-of-freedom  emerge  at  fine  scale  associated  with  our  correction  to  the 
Cauchy  Born  approximation,  v. 

This  key  difference  from  other  methods  holds  potential  for  further  development.  We 
speculate  that  difficulties  in  the  literature  for  establishing  matching  conditions  for  general 
thermodynamic  quantities  across  scales  can  perhaps  be  obviated  by  such  an  approach.  We 
contrast  this  to  the  displacement  matching  and  ghost  force  correction  of  the  local-nonlocal 
quasi-continuum  formulation  (Shenoy  et  al.,  1999),  and  the  lattice  summation  rules 
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invoked  in  the  nonlocal  quasi-continuum  method  (Knap  and  Ortiz,  2001),  both  of  which 
require  atomic  and  continuum  kinematics  to  be  described  by  the  same  coordinate  system. 
Likewise,  we  point  to  the  bridging  scales  method  (Park  et  ah,  2005)  that  employs  a  matrix 
containing  values  of  the  finite  element  shape  functions  evaluated  at  all  the  atomic  positions 
within  the  finite  element,  a  technique  again  requiring  that  the  scales  share  the  same 
coordinate  system. 

The  asymptotic  nature  of  our  method  also  permits  the  consideration  of  larger  scales  of 
the  continuum  relative  to  the  fine  scale  coordinates,  though  we  note  that  the  present  paper 
focuses  on  the  atomistic  influence  upon  the  response  of  a  single  continuum  material 
element  as  opposed  to  emphasizing  solution  of  the  macroscopic  equation  (53).  For 
example,  our  method  could  easily  be  incorporated  numerically  in  a  mixed-mesh  scenario  in 
which  only  regions  of  special  interest  are  addressed  using  multiscale  finite  elements 
incorporating  atomistic  homogenization,  with  the  remainder  of  the  domain  discretized 
with  standard  continuum  elastic  or  elastic-plastic  finite  elements.  Similar  techniques  have 
been  proposed  in  the  context  of  two-scale  continuum  homogenization  of  solid  mixtures 
(Oden  et  ah,  1999;  Ghosh  et  ah,  2001). 

Outside  of  the  context  of  the  homogenization  equations,  with  the  attendant  material 
embedding  and  convergence  features,  the  very  same  results  can  be  obtained  without  the 
present  developments.  Indeed,  our  intent  has  been  to  show  validation  of  the  method  with 
conjugate-gradient-based  lattice  statics  calculations  (part  (b)  of  Figs.  4-7),  and  then  to 
show  how  one  can  extend  the  calculations  to  the  setting  of  finite  deformation  plasticity. 
Although  the  kinetics  are  prescribed  a  priori,  a  limitation  discussed  above,  to  date  we  are 
unaware  of  another  atomistic-continuum  method  enabling  all  four  of  the  following 
features  in  the  context  of  multiscale  modeling  of  defective  crystals:  (i)  a  periodicity 
assertion  effortlessly  addressing  the  multiscale  representation  of  a  finite  defect  density  (i.e., 
non-isolated  defects);  (ii)  governing  equations  self-consistently  derived  over  two  length 
scales;  (hi)  an  explicit  mathematical  statement  for  the  corrected  atomic  displacements  near 
defects  parametrically  dependent  on  hyperelastic  lattice  distortions;  and  (iv)  an  evolving 
intermediate  configuration  accounting  for  the  stress-free  part  of  dislocation  motion  in  a 
general  sense,  independent  of  the  choice  of  kinetic  model  of  plasticity. 

8.  Conclusions 

A  multiscale  method  predicated  on  asymptotic  homogenization  as  been  developed  and 
implemented.  The  technique  enables  computation  of  the  macroscopic  (i.e.,  homogenized) 
mechanical  properties  such  as  effective  elastic  stiffness  and  net  stress  of  a  deforming  unit 
cell  consisting  of  periodically  arranged  discrete  atoms  at  the  microscale.  The  formulation 
directly  accounts  for  the  effects  of  defects  (e.g.  vacancies,  dislocations,  and/or 
misorientation  boundaries)  on  the  homogenized  mechanical  properties,  as  well  as  the 
lattice-preserving  kinematics  of  finite  plastic  deformation  associated  with  dislocation  flux, 
for  example,  all  as  functions  of  the  state  of  deformation.  From  a  computational 
standpoint,  our  homogenization  method  has  been  demonstrated  to  be  an  accurate 
alternative  to  incremental  conjugate  gradient  schemes  in  terms  of  prediction  of  minimum 
energy  configurations  of  atomic  degrees  of  freedom  in  statically  deforming  lattices 
containing  defects.  The  model  framework  enables  consideration  of  large  continuum  scale 
elements  embedded  with  finite  defect  densities,  the  latter  represented  explicitly  by 
periodically  perturbed  atomic  arrangements  at  the  fine  scale. 
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The  homogenization  method  was  implemented  numerically  (specifically  here  in  single- 
finite-element  integration  point  simulations)  to  study  the  nonlinear  elastic  response  of  BCC 
tungsten  lattices  containing  periodically  distributed  vacancies,  screw  dislocations  and 
dipoles,  and  low-angle  twist  boundaries.  It  was  found  that  strain  energies  associated  with 
vacancies,  screw  dislocations,  and  screw  dislocation  dipoles  tended  to  increase  with  applied 
uniaxial  stretching,  while  strain  energies  of  twist  boundaries  tended  to  decrease  with 
stretch.  Elastic  stiffness  in  the  direction  of  stretch  tended  to  decrease  with  increasing 
dislocation  content,  and  increase  with  twist  grain  boundary  area.  Anisotropy  of  the  elastic 
constants  of  tungsten,  nominally  isotropic,  was  also  demonstrated  at  applied  deformations 
and  in  the  presence  of  defects.  The  model  was  implemented  in  a  limited  fashion  to  study 
the  elastic-plastic  response  of  tungsten  lattices  containing  fixed  distributions  of  screw 
dislocations.  In  pure  shear  simulations,  influences  of  dislocation  density  on  shear  modulus 
and  energy  density  were  apparent. 

Finally,  we  summarize  the  limitations  of  our  approach.  Firstly,  due  to  the  periodicity 
requirement  imposed  by  the  homogenization  scheme,  only  perfectly  periodic  defect 
arrangements  may  be  modeled;  isolated  defects  in  infinite  lattices,  while  not  considered  at 
present,  could  conceivably  be  modeled  upon  substantial  modifications  of  the  theory  and 
implementation.  Secondly,  our  capability  is  currently  limited  to  addressing  purely 
mechanical,  rate-independent  responses,  with  the  microscopic  description  predicated  upon 
a  static  representation.  Hence,  finite  temperatures  and  viscous  behavior  are  not  captured, 
and  plasticity  kinetics  (e.g.,  thermally  assisted  dislocation  motion)  must  be  prescribed  a 
priori.  Defect  densities  and  numbers  of  atoms  remain  fixed  throughout  our  fine  scale 
simulations,  meaning  such  complex  phenomena  as  vacancy  migration,  dislocation 
generation,  and  dislocation  annihilation  are  not  considered  explicitly  in  the  atomistic 
domain.  Generalizations  to  the  model  framework  to  address  the  above  limitations  remain 
to  be  considered  in  future  work. 
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